Colorado  School  of  Mines 

Golden,  Colorado  80401 


Center  for  Wave  Phenomena 


Migration  Velocity  Analysis 


Zhenyue  Liu 


Colorado  School  of  Mines 

Golden,  Colorado  80401 


Center  for  Wave  Phenomena 


3rATafFriT~i~ 


^  puBiic  reieoMI 

. . .  Uaixmiwd 


) 


iMt  aJB?fieTEi>  a 


^997(I7U  134 


REPORT  DOCUMENTATION  PAGE 


0M8  No.  0704^0? 88 


'■foomog  Durden  ^or  mu  conecTton  of  information  l^  ntimateo  to  average  t  ngyr  resoonse.  inciuoina  the  time  for  reviewma  in^trurtinn*  • 

gathering  and  maintaining  me  data  needed,  and  comoieting  and  reviewing  the  coitection  of  information  Send  commenw  regarding  tnu  burden  estimate  or  anv^othM 
collection  of  information,  -nciuding  auggeationt  for  reducing  tm^  ouroen  to  Waihinoton  neadduanera  Services  Dirertorarefnr  intnrm*»^  other  ausect  of  thit 


1.  AGENCY  USE  ONLY  (Leave  biank)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  OATES  COVERED 

March  1995  Type  A:  I/95-3/95  Interim 


4.  TITLE  AND  SUBTITLE 

Migration  Velocity  Analysis 


6.  AUTHOR(S)  — - 

Zhenyue  Liu 


7.  PERPORMING  ORGANIZATION  NAME(S)  AND  AOOR£SS(ES) 

Center  for  Wave  Phenoinena 
Colorado  School  of  Mines 
Golden,  CO  80401 


9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  AODRESS(ES) 

Mathematical,  Computer  &  Information  Sciences 
Division;  Applied  Analysis  Program 


5.  PUNOING  NUMBERS 

Grant  #N00014-9 1-J- 
1267 

R&T  Project  # 
1112669 - 08 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


CWP-168 


10.  SPONSORING/ MONITORING 
AGENCY  REPORT  NUMBER 


11.  SUPPLEMENTARY  NOTES  . .  '  . 

Ph.D.  thesis.  Mathematical  &  Computer  Sciences,  Colorado  School  of 
Mines 

12i.  OISTRIBUTIQM/ AVAILABILITY  ^TATCmpmt  _ , _ _ I  Uk  niETniBUTION  rnnr 


unclassified 


f  pmfiJ2lTiOE  I 

j  &B|:>zov€d  tcu  .gatiiic  ! 

Pise  Ibunciea; 


13.  ABSTRACT  (Maximum  200  words) 

This  thesis  addresses  the  development  of  a  general,  quantitative  representation  between  the 
imaged  depths  and  migration  velocity.  Based  on  ray  theory,  I  have  found  that  the  migration  process 
can  be  described  by  using  the  so-called  imaging  equations.  From  this  viewpoint,  I  have  studied 
relationships  between  the  residual  moveout  and  the  error  of  migration  velocity  in  a  general  context. 
To  handle  lateral  velocity  variations,  I  used  a  perturbation  method  to  derive  an  analytical  formula  to 
iteratively  update  the  velocity  from  residual  moveout.  My  formula  estimates  the  update  in  velocity 
by  computing  a  derivative  function  of  imaged  depths  with  respect  to  velocity  in  a  general  background 
medium  context.  This  formula  is  more  accurate  than  conventional  formulas  based  on  hyperbolic 
residual  moveout  when  the  medium  has  strongly  lateral  velocity  variations.  In  addition,  I  discuss 
the  sensitivity  of  migration-based  velocity  analysis,  and  show  what  factors  affect  the  sensitivity. 

The  theory  and  formulas  developed  here  are  available  both  for  the  2-D  case  and  for  the  3-D  case. 
The  methodology  in  this  thesis  have  been  tested  on  synthetic  data  (including  the  Marmousi  data) 
and  on  physical-tank  data  that  are  recorded  in  a  scaled  earth  model. 


15.  NUMBER  OF  PAGES 

88 


16.  PRICE  CODE 


14.  SUBJECT  TERMS 

migration 

;  velocity  analysis 

;  imaging  equations 

17.  SECURITY  CLASSIFICATION 

OF  REPORT 

Unclassified 

18.  SECURITY  CLASSIFICATION 

OF  THIS  PAGE 

Unclassified 

19.  SECURITY  CLASSIFICATION 

OF  ABSTRACT 

Unclassified 

NSN  7540-0! -280-5500 


DTIC  QUALITY  mSPECTED  3 


Standard  Form  298  (Rev  2*89) 

Pr^nbed  by  ANSI  Std  Z39*i8 
298-102 


DEPARTMENT  OF  THE  NAVY 

OFFICE  OF  NAVAL  RESEARCH 
SEATTLE  REGIONAL  OFFICE 
1107  NE  45TH  STREET.  SUITE  350 

SEATTLE  WA  98105-4631  IN  REPLY  REFER  TO: 


4330 
ONR  247 
1 1  Jul  97 


From:  Director,  Office  of  Naval  Research,  Seattle  Regional  Office,  1 107  NE  45th  St.,  Suite  350, 
Seattle,  WA  98105 

To:  Defense  Technical  Center,  Attn:  P.  Mawby,  8725  John  J.  Kingman  Rd.,  Suite  0944, 

Ft.  Belvoir,VA  22060-6218 

Subj:  RETURNED  GRANTEE/CONTRACTOR  TECFINICAL  REPORTS 


1.  This  confirms  our  conversations  of  27  Feb  97  and  1 1  Jul  97.  Enclosed  are  a  number  of 
technical  reports  which  were  returned  to  our  agency  for  lack  of  clear  distribution  availability 
statement.  This  confirms  that  all  reports  are  unclassified  and  are  “APPROVED  FOR  PUBLIC 
RELEASE”  with  no  restrictions. 


2.  Please  contact  me  if  you  require  additional  information.  My  e-mail  is  silverr@onr.navy.mil 
and  my  phone  is  (206)  625-3196. 

/"■  -A/. 

ROBERT  J.  SILVERMAN 


CWP-168 
March  1995 


Migration  Velocity  Analysis 

Zhenyue  Liu 

—  Doctoral  Thesis  — 

Mathematical  and  Computer  Sciences 


Center  for  Wave  Phenomena 
Colorado  School  of  Mines 
Golden,  Colorado  80401 
303/273-3557 


ABSTRACT 


Prestack  depth  migration,  which  can  both  handle  dipping  reflectors  and  lateral  ve¬ 
locity  variations,  is  a  robust  method  for  imaging  complex  structures.  In  order  to  process 
data  by  this  method,  one  often  needs  to  have  a  more  accurate  velocity  model  than  may 
be  obtained  from  simple  velocity  analysis  methods,  such  as  normal  moveout. 

On  the  other  hand,  prestack  depth  migration  itself  is  an  attractive  tool  for  doing 
velocity  analysis  because  of  its  high  sensitivity  to  the  velocity  model.  Two  approaches 
to  migration  velocity  analysis  have  become  prominent:  depth-focusing  analysis  (DFA) 
and  residual-curvature  analysis  (RCA).  So  far,  the  formulas  used  in  both  approaches 
for  updating  velocity  are  derived  under  the  assumptions  of  horizontal  reflectors,  lateral 
velocity  homogeneity,  or  small  offset.  Therefore,  those  conventional  approaches  lack 
sufficient  computational  efficiency  and  accuracy  when  a  velocity  distribution  has  strong, 
lateral  variations. 

This  thesis  addresses  the  development  of  a  general,  quantitative  representation  be¬ 
tween  the  imaged  depths  and  migration  velocity.  Based  on  ray  theory,  I  have  found 
that  the  migration  process  can  be  described  by  using  the  so-called  imaging  equations. 
From  this  viewpoint,  I  have  studied  relationships  between  the  residual  moveout  and  the 
error  of  migration  velocity  in  a  general  context.  For  a  medium  with  weak  lateral  veloc¬ 
ity  variations,  I  have  proved  that  the  residual  velocity  can  be  estimated  by  hyperbolic 
residual  moveout  and  no  iteration  is  required.  To  handle  lateral  velocity  variations,  1 
used  a  perturbation  method  to  derive  an  analytical  formula  to  iteratively  update  the 
velocity  from  residual  moveout.  In  my  derivation,  I  impose  no  limitation  on  offset,  dip, 
or  velocity  distribution  as  long  as  the  velocity  perturbation  is  sufficiently  small.  My 
formula  estimates  the  update  in  velocity  by  computing  a  derivative  function  of  imaged 
depths  with  respect  to  velocity  in  a  general  background  medium  context.  This  formula 
is  more  accurate  than  conventional  formulas  based  on  hyperbolic  residual  moveout  when 
the  medium  has  strongly  lateral  velocity  variations;  therefore,  it  is  a  major  contribution 
of  this  thesis.  Based  on  this  formula,  I  revise  the  RCA  method  for  doing  velocity  analysis 
in  complex  structural  media.  In  addition,  with  help  of  the  imaging  equations,  I  analyse 
properties  in  a  common  image  gather  (CIG),  discuss  the  sensitivity  of  migration-based 
velocity  analysis,  and  show  what  factors  affect  the  sensitivity. 

The  theory  and  formulas  developed  here  are  available  both  for  the  2-D  case  and  for 
the  3-D  case.  For  converted  waves  and  anisotropic  media,  I  also  offer  some  suggestions 
on  how  to  develop  corresponding  formulas. 

The  theory  and  methodology  in  this  thesis  have  been  tested  on  synthetic  data  (in¬ 
cluding  the  Marmousi  data)  and  on  physical-tank  data  that  are  recorded  in  a  scaled  earth 
model. 
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Chapter  1 


INTRODUCTION 


Structural  imaging  is  an  essential  process  in  seismic  exploration.  A  successful  imple¬ 
mentation  of  this  process  relies  on  knowledge  of  background  parameters  for  propagation 
of  data  into  the  earth.  Among  these  parameters,  velocity  is  a  critical  one  because  of  its 
coherent  relation  with  the  dynamic  characteristics  of  reflection  signals.  Therefore,  ve¬ 
locity  estimation  has  become  one  of  the  most  important  tasks  in  exploration  seismology. 
Velocity  analysis  is  the  process  of  estimating  background  velocities  from  measurements  of 
moveout  or  stacking  power.  In  general,  there  are  three  ways  for  doing  velocity  analysis: 
normal  moveout,  dip  moveout  and  prestack  migration. 

Normal  moveout  (NMO)  is  the  most  basic  method  for  determining  velocities  from 
seismic  data.  However,  this  method  assumes  horizontal  reflectors  and  lateral  velocity 
homogeneity.  Dip-moveout  (DM0)  improves  on  NMO  processing  by  properly  accom¬ 
modating  dipping  reflectors,  but  still  cannot  accommodate  lateral  velocity  variations. 
Compared  to  these  two  methods,  prestack  migration  provides  a  powerful  tool  for  doing 
velocity  analysis  because  of  its  high  sensitivity  to  the  velocity  error  and  its  ability  to 
handle  both  reflector  dips  and  lateral  velocity  variations. 

1.1  Paper  Review 

Two  approaches  to  migration  velocity  analysis  have  been  developed:  depth-focusing 
analysis  (DFA)  and  residual-curvature  analysis  (RCA).  Depth-focusing  analysis  is  based 
on  using  stacking  power  to  measure  velocity  error.  Residual-curvature  analysis  is  based 
on  using  residual  moveout  to  measure  velocity  error.  Here,  I  will  give  an  overview  on 
papers  of  these  two  approaches.  Moreover,  I  will  make  a  brief  comparison  between 
migration-based  analysis  and  the  reflection  tomography  approach. 

1.1.1  Depth- focusing  analysis 

The  idea  of  velocity  analysis  based  on  dififerential  solutions  of  the  scalar  wave  equa¬ 
tion  was  first  introduced  by  Doherty  and  Claerbout  (1974).  They  used  the  15-degree 
finite-difference  migration  algorithm  and  worked  with  single  CMP  gathers  to  carry  out 
velocity  estimation  from  unstacked  seismic  data.  Yilmaz  and  Chambers  (1984)  extended 
the  migration-based  approach  to  include  more  than  one  CMP  gather  in  each  analysis. 
This  extension  allows  proper  treatment  of  dipping  events  and  yields  velocity  information 
that  is  more  appropriate  for  use  in  migration.  The  migration  process  was  implemented 
by  using  downward  extrapolation  in  the  Fourier  transform  domain.  Jeanot  et  al.  (1986) 
then  generalized  the  method  to  prestack  depth  migration  in  the  form  of  depth-focusing 
analysis  by  determining  the  focused  depth  in  the  process  of  downward  continuation.  Sub- 
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sequently,  MacKay  and  Abma  (1989)  showed  how  the  DFA  data  volume  could  be  used  to 
form  a  well-focused  seismic  section.  The  new  section,  the  focal-surface  image,  estimates 
the  results  of  prestack  depth  migration  using  the  updated  velocities. 

The  DFA  approach  is  based  on  these  principles:  (1)  when  the  migration  velocities 
are  exact,  the  two  imaging  conditions,  zero  time  and  zero  offset,  yield  a  focused  image 
during  downward  continuation;  (2)  when  the  migration  velocities  are  in  error,  reflected 
energy  collapses  to  zero  offset  at  depths  that  are  inconsistent  with  the  zero-time  imaging 
condition;  (3)  by  interpreting  the  nonzero  times  at  which  focusing  actually  occurs,  the 
migration  velocities  can  be  updated  iteratively. 

A  simple  formula  is  used  to  update  velocity  in  the  DFA  method: 

K  =  (1.1) 

where  is  the  migration  depth,  Kn  is  the  migration  velocity,  6  is  half  the  depth- 
difference  between  the  focusing  and  migration  depth,  and  is  the  updated  velocity.  The 
formula  (1.1)  is  strictly  valid  for  horizontal  reflectors,  small  offset  angles,  and  constant 
velocities.  For  a  depth-variable  velocity,  Vr  and  Vm  are  viewed  as  the  rms  velocities.  If 
the  effect  of  steep  dip  is  included,  the  modified  formula  (MacKay,  1991)  is 


where  6m  is  the  migrated  dip  and  Or  is  the  real  structural  dip.  Unfortunately,  this  formula 
cannot  be  used  to  update  velocities  because  9r  is  not  known. 

MacKay  and  Abma  (1992)  discussed  how  to  avoid  divergence  in  the  velocity-updating 
procedure,  when  reflector  dips  and  lateral  velocity  variations  exist.  They  applied  a  damp>- 
ing  factor  to  the  interpreted  depth  errors.  In  addition,  the  interpretation  of  the  focusing 
depth  is  frequently  affected  by  seismic  energy  from  dipping  interfaces,  diffractions,  and 
noise.  To  reduce  possible  ambiguity  in  the  DFA  interpretation  process,  MacKay  and 
Abma  (1993)  introduced  a  new  attribute  for  determining  focusing  based  on  estimates  of 
the  radius  of  wavefront  curvature. 

The  large  amount  of  the  ion  required  for  downward  continuation  may  limit  use  of  the 
DFA  approach  for  routine  processing.  Kim  and  Gonzalez  (1991)  presented  the  Kirchhoff 
integral  approach  to  the  downward  continuation.  The  flexibility  of  the  Kirchhoff  integral 
allows  one  to  compute  only  the  zero-offset  trace  at  each  depth  point  and  to  avoid  most 
of  computation  for  the  downward  continuation  of  unstacked  data. 


1.1.2  Residual-curvature  analysis 

The  formula  for  updating  velocity  in  the  DFA  approach  is  dip-limited,  so  repeated 
prestack  migration  is  required  even  for  a  constant-velocity  medium.  In  comparison, 
residual-curvature  analysis  (RCA)  is  an  alternative  to  migration  velocity  analysis  that  is 
able  to  overcome  the  dip  limitations  of  DFA. 

In  the  RCA  method,  the  migrated  prestack  data  are  sorted  into  common-image  gath- 
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ers  (CIGs).  They  have  also  been  called  common-receiver  gathers  (CRGs)  or  common- 
reflection  point  gathers  (CRPs).  In  each  common-image  gather,  the  migrated  data  have 
the  same  imaged  horizontal  location.  Residual-curvature  analysis  is  based  on  this  princi¬ 
ple;  after  prestack  migration  with  the  correct  velocity,  the  imaged  depths  in  a  GIG  must 
be  the  same  regardless  of  offset;  otherwise,  after  prestack  migration  with  an  erroneous 
velocity,  the  imaged  depths  in  a  CIG  from  different  offsets  will  differ  from  each  other; 
the  differences  of  imaged  depth  in  CIGs  provide  information  for  updating  the  velocity 
iteratively. 

Al-Yahya  (1989)  discussed  residual-curvature  analysis  by  iterative  shot  profile  mi¬ 
gration.  He  measured  the  velocity  errors  by  estimating  the  curvature  of  residual  moveout. 
The  residual  moveout  in  a  CIG  in  his  paper  is  represented  by 

4  =  +  (7"  -  l)«^  (1.3) 

where  a  is  the  offset  between  the  shot  and  the  imaged  location,  is  the  imaged  depth, 
z  is  the  actual  depth,  and  7  is  the  ratio  of  the  true  rms  slowness  to  the  migration  rms 
slowness.  The  derivation  of  equation  (1.3)  is  based  on  small  offset,  horizontal  reflectors, 
and  laterally  invariant  velocity.  Later,  Lee  and  Zhang  (1992)  generalized  the  condition 
of  horizontal  reflector  into  small-angle  dip  by  using  the  dip-corrected  residual  moveout 
equation,  from  which  equation  (1.3)  is  modifled  by 

=  7^2^  -I-  (72  -  l)(a2  -f  (a  -  Uo)^),  (1.4) 

where  uq  is  a  quantity  related  to  the  reflector  dip.  With  the  help  of  this  more  accurate 
equation,  one  can  implement  residual  migration  to  avoid  remigration,  or  reduce  the 
number  of  the  number  of  iterations  required  to  produce  a  satisfactory  image.  However, 
Lee  and  Zhang’s  approach  involves  two  unknown  parameters — the  true  velocity  and  the 
reflector  dip,  at  each  imaged  location,  which  makes  this  approach  less  convenient  than 
the  determination  of  one  single  parameter. 

Residual  moveout  in  common  offset  migration  is  relatively  independent  of  reflector 
dip,  so  common  offset  migration  provides  a  better  approach  in  residual-velocity  analysis 
than  common  shot  migration.  Deregowski  (1990)  realized  this  advantage  of  common  off¬ 
set  migration  over  other  migration  techniques  and  developed  a  velocity  analysis  method 
by  using  common-offset  depth  migration.  In  his  approach,  the  velocity  is  updated  it¬ 
eratively  by  using  the  residual-moveout  correction.  Liu  and  Bleistein  (1992)  derived  a 
general  residual-moveout  representation  under  the  assumption  of  small  offset  and  proved 
that  residual  moveout  in  common  offset  migration  is  insensitive  to  reflector  dip  for  a 
depth-dependent  velocity.  This  conclusion  implies  the  velocity  can  be  corrected  without 
a  further  iteration  for  a  laterally  invariant  velocity.  To  handle  complex  structures,  Lan- 
ford  and  Levander  (1993)  proposed  an  RCA  method  that  is  based  on  a  layer-stripping 
Kirchhoff  algorithm  with  geometrical  ray  tracing  in  heterogeneous  media.  They  assumed 
that  the  medium  is  composed  of  constant- velocity  layers  and  updated  velocity  along  the 
ray  path  in  a  way  similar  to  tomographic  reconstruction.  By  using  perturbation  theory, 
Liu  and  Bleistein  (1994)  developed  an  analytical  formulation  for  residual  moveout  that 
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is  valid  for  any  velocity  distribution  and  any  subsurface  structure.  This  formulation  can 
also  be  used  to  analyze  the  sensitivity  of  migration  velocity  analysis. 

1.1.3  Comparison  with  reflection  tomography 

As  an  alternative  to  migration  velocity  analysis,  reflection  tomography  is  an  impor¬ 
tant  approach  to  velocity  analysis  that  handles  lateral  velocity  variations  for  velocity 
estimation.  Therefore,  it  is  essential  to  make  a  comparison  between  migration  velocity 
analysis  and  reflection  tomography. 

The  tomography  technique  (refraction  tomography)  has  been  widely  used  in  Bore¬ 
hole  Geophysics  and  Vertical  Seismic  Profile  (VSP)  to  estimate  medium  velocity  based 
on  information  of  the  first  arrival  traveltime.  It  is  called  reflection  tomography  when  this 
technique  is  applied  to  surface  seismic  reflection  data.  In  the  tomographic  approach  to 
the  inverse  problem,  the  medium  velocity  is  determined  by  minimizing  a  misfit  function 
that  represents  the  deviation  of  traveltimes.  The  formulation  in  the  tomographic  ap¬ 
proach  is  set  up  analytically,  so  one  can  analyse  some  aspects  of  velocity  inversion,  such 
as  stability,  resolution  and  computational  efficiency. 

However,  there  are  some  difficulties  in  the  present  reflection  tomographic  approach. 
First,  picking  reflected  events  from  seismic  signals  in  reflection  tomography  is  much  more 
difficult  than  picking  the  first  arrivals  in  refraction  tomography.  The  events  in  the  unmi¬ 
grated  data  are  frequently  distorted  by  correlated  noise.  This  problem  can  be  improved 
by  picking  reflectors  in  the  postmigrated  domain  (Stork,  1992).  Secondly,  the  inversion 
problem  in  reflection  tomography  is  very  ill-posed.  Compared  to  refraction  tomography, 
the  raypath  coverage  in  reflection  tomography  is  more  sparse  and,  therefore,  uniqueness  is 
an  issue,  with  multiple  solutions  being  a  possibility.  Finally,  the  construction  of  raypaths 
(actually,  specular  raypaths)  strongly  relies  on  the  reflector  position,  both  on  depths  and 
on  slopes.  Although  the  depths  can  be  included  in  the  inverse  problem  as  a  parameter  to 
be  determined  (Stork  and  Clayton,  1991),  generally  it  is  difficult  to  determine  the  slopes 
accurately. 

In  contrast,  migration  velocity  analysis  can  use  the  semblance  method  to  measure 
residual  moveouts  and  pick  reflectors  in  the  stacked,  migrated  domain.  In  this  way,  mi¬ 
gration  velocity  analysis  will  give  more  reliable  information  for  picking.  Furthermore, 
with  the  help  of  the  macro  model  concept,  migration  velocity  analysis  can  apply  a  re¬ 
cursive  procedure  to  estimate  interval  velocities.  Thus,  the  computation  can  be  made 
more  stable  by  applying  the  method  to  only  a  few  unknowns  at  each  recursion.  Finally, 
formulation  in  migration  velocity  analysis  can  be  set  up  analytically  by  using  pertur¬ 
bation  theory  which  is  similar  to  the  analysis  done  in  tomography  (Liu  and  Bleistein, 
1994).  The  estimate  of  velocities  in  this  formulation  also  depends  on  specular  raypaths, 
as  in  reflection  tomography.  In  contrast,  Liu  and  Bleistein’s  approach  selects  the  specular 
raypaths,  based  on  the  stationary-phase  principle  that  produces  the  migration  imaging, 
instead  of  on  the  measurement  of  reflector  slopes.  Therefore,  this  approach  is  more  stable 
and  less  computationally  costly  than  reflection  tomography. 

In  summary,  both  migration  velocity  analysis  and  reflection  tomography  are  impor¬ 
tant  methods  for  handling  lateral  velocity  variation  in  velocity  estimation.  The  similarity 
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between  these  two  methods  allows  them  to  “learn”  from  each  other  in  order  to  solve  the 
problem  of  velocity  estimation  in  complex  structures.  For  instance,  Liu  and  Bleistein 
(1994)  derived  an  analytical  formula  to  describe  the  relation  between  residual  moveout 
and  residual  velocity,  the  idea  for  which  was  suggested  by  tomography.  On  the  other 
hand,  reflection  tomography  may  improve  stability  by  using  the  macro  model  concept. 

1.2  Limitations  in  Conventional  Approaches 

Migration  velocity  analysis  must  address  two  issues  to  succeed:  (1)  how  to  establish 
a  criterion  for  knowing  if  a  migration  velocity  is  acceptable;  (2)  how  to  update  the  ve¬ 
locity,  if  it  is  unacceptable.  In  DFA,  a  migration  velocity  is  acceptable  if  the  difference 
between  migration  depth  and  focusing  depth  is  zero.  In  RCA,  a  migration  velocity  is 
acceptable  if  the  difference  between  imaged  depths  from  different  offsets  is  zero.  Quanti¬ 
tatively,  these  differences  can  be  used  to  update  the  velocity  in  a  routine  way.  To  date,  a 
variety  of  updating  formulas  have  been  developed.  These  formulas  degrade  with  increas¬ 
ing  complexity  of  media  (lateral  velocity  variation  and  reflector  dip),  because  of  rough 
approximations  used  to  estimate  velocity.  Although  iteration  generally  is  helpful  in  ob¬ 
taining  a  more  accurate  velocity,  a  good  approximation  not  only  reduces  the  number  of 
iterations  but  also  guarantees  the  convergence.  The  conventional  formulas  for  updating 
velocity  were  derived  under  one  or  more  assumptions  as  follows: 

(i)  lateral  velocity  homogeneity; 

(ii)  small  offset; 

(iii)  horizontal  reflector. 

In  MacKay  and  Abma’s  DFA  approach  (1992),  all  three  assumptions  are  used.  In  RCA, 
the  assumptions  vary  with  migration  types.  For  common-shot  migration  or  common- 
receiver  migration,  Al-Yahya’s  formula  (1989)  used  all  three  assumptions.  Lee  and  Zhang 
(1992)  derived  a  formula  that  replaces  assumption  (iii)  by  the  assumption  of  small  dip. 
For  common-offset  migration,  Deregowski’s  formula  (1990)  used  the  first  two  assump¬ 
tions.  [This  statement  was  verified  by  Liu  and  Bleistein  (1992).]  Both  DFA  and  RCA 
use  assumptions  (i)  and  (ii).  Under  the  assumption  of  a  small  offset,  the  residual  moveout 
can  be  approximated  by  a  hyperbola  or  a  parabola.  However,  this  approximation  is  poor 
when  the  velocity  has  an  obvious  lateral  variation. 

Conventional  approaches  for  updating  the  velocity  inspect  information  at  each  mid¬ 
point  individually.  This  information  is  not  reliable  when  conflicting  events  exist  because 
velocity  estimation  will  depend  on  which  event  is  used.  Furthermore,  the  conventional 
approaches  update  the  average  velocities  that  are  assumed  to  be  equal  to  the  rms  veloc¬ 
ities  and  then  convert  rms  velocities  into  interval  velocities  by  the  Dix  equation  or  other 
methods.  Unfortunately,  when  the  velocity  has  a  lateral  anomaly,  this  kind  of  average 
velocity  may  differ  significantly  from  the  rms  velocity  (Lynn  and  Claerbout,  1982;  Toldi, 
1989),  which  may  cause  divergence  in  the  velocity-updating  procedure. 

Because  of  such  limitations,  conventional  approaches  to  velocity  analysis  encounter 
difficulties  in  handling  complex  structures  (such  as  the  Marmousi  model).  To  overcome 
these  limitations,  some  geophysicists  used  a  macro-model  (layer-stripping)  method  for 
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estimating  the  interval  velocity  directly.  A  macro-model  consists  of  velocities  and  veloc¬ 
ity  interfaces  (subsurfaces).  The  model  is  determined  by  layer-stripping,  in  a  top  down 
procedure;  the  velocity  of  each  layer  is  updated  iteratively  by  using  depth-focusing  anal¬ 
ysis  or  residual-curvature  analysis;  the  velocity  interface  is  imaged  by  using  the  corrected 
velocity.  The  whole  process  may  be  accomplished  with  interpretation  based  on  geological 
knowledge.  So  far,  the  present  macro-model  method  has  two  significant  drawbacks.  First, 
the  velocity  distribution  is  assumed  to  be  constant  in  layers.  In  fact,  velocity  may  vary 
within  a  layer.  Secondly,  there  has  not  been  an  efficient  formula  for  updating  velocity, 
which  results  in  inaccuracy  and  inefficiency  in  velocity  estimates. 

In  summary,  prestack  depth  migration  is  sensitive  to  the  velocity  model.  Residual 
curvatures  in  migrated  data  provide  information  indicating  if  the  migration  velocity  is 
correct  or  not.  However,  there  exist  drawbacks  in  the  present  formulas  for  updating 
velocity  due  to  rough  approximations  used.  Therefore,  the  derivation  of  a  more  accurate 
formula  is  required  to  handle  complex  structures. 

1.3  Technical  Work  in  this  Thesis 

Migration  imaging  depends  on  what  velocity  is  used  and,  on  the  other  hand,  migra¬ 
tion  imaging  determines  what  velocity  should  be  used.  Therefore,  it  is  essential  to  give 
a  quantitative  description  to  the  relationship  between  migration  imaging  and  migration 
velocity,  especially  to  the  relationship  between  residual  moveout  and  residual  velocity. 

In  Chapter  2,  I  set  up  the  imaging  equations  which  are  a  kinematic  representation  of 
migration  processes.  I  address  the  equivalency  of  the  imaging  equations  to  Snell’s  Law, 
so  these  equations  also  can  be  applied  to  converted  waves  and  anisotropic  media.  The 
imaging  equations  are  the  foundation  of  the  velocity  analysis  formulas  derived  in  Chapter 
4  and  Chapter  5. 

In  Chapter  3,  I  analyze  image  properties  at  a  common  image  gather  based  on  the 
imaging  equations.  For  constant  layered  velocity,  I  derive  formulas  for  the  sensitivity  of 
velocity  analysis  and  show  what  factors  affect  the  accuracy  of  migration  velocity  analysis. 

In  Chapter  4,  I  derive  the  residual  moveout  representations  under  the  assumption 
of  small  offset.  Based  on  these  representations,  I  define  the  residual  moveout  (RMO) 
velocity  and  show  how  the  RMO  velocity  is  related  to  the  rms  velocity.  I  prove  that  the 
RMO  velocity  can  be  estimated  by  hyperbolic  residual  moveout  without  iteration  and  it 
is  a  good  approximation  of  the  rms  velocity  for  a  medium  with  weakly  lateral  velocity 
variations.  In  addition,  I  show  a  failure  example  on  a  simple  two-layer  model  where  the 
lateral  variation  is  no  longer  small. 

In  Chapter  5,  I  develop  iterative  approaches  to  handle  lateral  velocity  variations. 
Based  on  the  perturbation  method,  I  derive  an  analytical  relationship  between  the  resid¬ 
ual  moveout  and  residual  velocity  and  discuss  some  computational  techniques  for  esti¬ 
mating  velocity  by  using  this  relationship.  My  formula  estimates  the  update  in  velocity 
by  computing  a  derivative  function  of  imaged  depths  with  respect  to  velocity  in  a  general 
background  medium  context.  This  formula  is  more  accurate  than  conventional  ones  based 
on  hyperbolic  residual  moveout  when  the  medium  has  strongly  lateral  velocity  variations. 
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In  addition,  I  propose  a  simple-iteration  algorithm  as  a  supplement  to  the  perturbation 
approach.  Furthermore,  I  derive  more  formulas  for  the  sensitivity  of  migration  velocity 
analysis,  which  are  a  development  of  the  formulas  in  Chapter  3. 

In  Chapter  6,  I  show  the  application  of  the  iterative  method  in  Chapter  5.  The 
computer  implementations  include  simple  synthetic  data,  physical-tank  data,  and  the 
Marmousi  data. 

The  formulas  of  velocity  analysis  in  this  thesis  are  derived  for  2-D  acoustic  media  as 
those  in  the  papers  Liu  and  Bleistein  (1992),  Liu  and  Bleistein  (1994),  Liu  and  Bleistein 
(1995).  Here,  I  have  partially  extended  these  results  to  3-D  and  I  also  present  suggestions 
on  how  to  extend  to  converted  waves  and  anisotropic  media. 
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Chapter  2 


IMAGING  EQUATIONS 


Let  us  assume  that  seismic  waves  are  generated  and  recorded  in  response  to  a  near¬ 
surface  source.  The  recorded  signals  will  provide  information  about  the  underground- 
velocity  distribution  and  reflector  geometry.  That  information  is  contained  in  amplitude 
and  phase  of  seismic  waves.  Under  the  assumption  of  high  frequency,  phase  information 
is  simplifled  to  travel  time  (the  WKBJ  approximation).  Compared  to  amplitude,  travel¬ 
time  information  in  seismic  waves  is  more  reliable.  For  example,  the  amplitude  in  a  P-P 
reflection  is  greatly  aflfected  by  the  existence  of  shear  waves,  but  the  traveltime  is  much 
less  affected.  Consequently,  studying  methods  based  on  traveltime  is  essential  in  explo¬ 
ration  seismology.  The  theory  that  deals  with  wave  propagation  for  very  high-frequency 
signals  is  called  geometrical  optics.  Using  this  tool,  Hubral  and  Krey  (1980)  did  funda¬ 
mental  work  for  measuring  traveltime  and  estimating  the  stacking  velocity.  In  this  thesis, 
I  am  also  going  to  utilize  this  tool  for  measuring  residual  moveout  and  estimating  the 
migration  velocity. 

Given  a  reflector  and  a  medium  velocity,  one  can  obtain  a  space-time  curve  in  the 
recorded  signals  by  ray  theory.  This  process  is  forward  modeling.  Conversely,  given  a 
time-space  curve  and  a  velocity,  one  can  determine  a  reflector  that  generates  this  space- 
time  curve.  This  process  is  migration  imaging.  Furthermore,  if  there  are  multichannel 
seismic  data,  velocity  itself  also  can  be  determined  from  the  relationship  between  the 
subsurface  imaging  and  the  migration  velocity.  This  process  is  velocity  analysis. 

In  this  chapter,  equations  for  describing  migration  imaging  will  be  set  up  based 
on  Snell’s  Law.  These  equations  are  the  kinematical  representations  of  the  migration 
processes. 

2.1  Mathematical  Derivation 

Consider  the  two  dimensional  situation.  I  denote  by  x  a  2-D  vector,  x  =  {x,z). 
Let  Xa  =  (rr,(^),Za(^))  be  source  positions  and  Xr  =  {xr{0,  Zr{^))  be  receiver  positions 
located  on  the  datum  surface  (not  necessarily  flat),  where  ^  is  a  position  parameter  on  the 
surface.  For  example,  in  common  offset,  ^  would  be  the  midpoint  of  a  source/receiver 
pair;  in  common  shot,  ^  would  be  a  coordinate  labeling  the  receiver  location,  that  is, 
representing  the  offset  from  the  fixed  source.  For  any  point  below  the  surface,  r,(Xj,x) 
denotes  traveltime  for  a  downgoing  wave  from  x^  to  x,  and  Tp(x,  x^),  denotes  traveltime 
for  an  upgoing  wave  from  x  to  x^.  (See  Figure  2.1.) 

Given  a  reflector,  one  can  compute  the  total  reflection  traveltime  t{^)  by  using 

r,(x„x)-l-rr(x,x,.)  =t(0,  (2.1) 
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datum  surface 


Fig.  2.1.  Reflection  geometry  for  the  2-D  case. 


where  the  reflection  point  x  satisfies  Snell’s  Law, 


[Vx  T^{Xs,x)  +  Vx  Tr{x,  X,.)]  •  —  =  0.  (2.2) 

Equations  (2.1)  and  (2.2)  describe  the  forward  modeling  process  for  a  given  reflec¬ 
tor,  so  they  are  called  the  modeling  equations.  On  the  other  hand,  the  total  reflection 
traveltime  function  t{^)  also  can  be  used  to  determine  the  reflector  position.  Notice 
that  the  reflection  point  x  =  {x,z)  is  dependent  on  the  position  ^  on  the  surface.  By 
differentiating  equation  (2.1)  with  respect  to  I  obtain 


dTg  dTj.  dx  dt 

gj- +  gj- +  [Vx  (T,  +  r.)J  ■  —  = 


(2.3) 


Using  equation  (2.2)  to  eliminate  the  third  term  on  the  left  side  of  equation  (2.3),  one 
obtains 


dT^{Xs,x)  dTr{x,Xr)  dt 


(2.4) 


Thus,  each  point  x  =  (x,  z)  of  the  reflector  (therefore,  the  reflector  itself)  is  determined 
by  the  simultaneous  solution  of  equations  (2.1)  and  (2.4);  i.e.,  x  will  satisfy 


T,{Xs,x)+Tr{x,Xr)  =t(0, 
dTg{Xa,x)  dTr{x,Xr)  _  dt 

ae  ae  ~  d(- 

Geometrically,  equation  (2.6)  shows  that  the  reflector  is  just  the  envelope  of  the  wave- 
fronts  from  source/receiver  pairs.  Therefore,  I  will  call  equation  (2.6)  the  envelope  con¬ 
dition.  This  condition  also  can  be  derived  by  using  the  stationary-phase  method  in  the 
asymptotic  expansion  of  the  Kirchhoff  migration  operator  to  reflection  data  under  the 
assumption  of  high  frequency.  However,  the  derivation  here  is  more  direct. 

Equations  (2.5)  and  (2.6)  are  the  dual  of  equations  (2.1)  and  (2.2).  They  describe 
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(2.5) 

(2.6) 


the  migration  imaging  process  for  a  given  space-time  curve,  so  I  call  them  the  imag¬ 
ing  equations.  These  two  equations  hold  for  any  kind  of  seismic  acquisition  geometry. 
Specifically,  let  us  consider  the  common  offset  case  with  a  fiat  datum  surface. 

With  no  loss  of  generality.  Let  Zg  ■=  =  0.  The  source  and  receiver  positions  can 

be  described  by  the  horizontal  coordinates  and  Xr  only.  Let  y  denote  the  midpoint 
and  h  the  half-offset: 

Xg  =  y-h,  Xr  =  y-\-h.  (2.7) 

Then,  since  ^  =  y,  equations  (2.5)  and  (2.6)  become 


T^{x„x)-{-Tr{x,Xr)  =t{y,h), 


(2.8) 


^  +  ^  =  -  f29) 

dy  ^  dy  ay- 

Equations  (2.8)  and  (2.9)  show  that  if  the  offset-traveltime  curve  t{y,  h)  (therefore,  dtfdy) 
is  known  for  given  h,  one  can  compute  an  imaged  depth  2  at  each  location  x.  Further¬ 
more,  if  a  migration  velocity  equals  the  true  velocity,  then  the  imaged  depth  2  must  be 
independent  of  the  offset  h;  otherwise,  for  a  wrong  migration  velocity,  2  will  vary  with 
the  offset  h.  Consequently,  the  imaged  depths  of  different  offsets  provide  information  for 
determining  the  velocity,  i.e.,  for  velocity  analysis. 

Similarly,  in  the  common-shot  case,  equations  (2.5)  and  (2.6)  become 


rs{Xs,x)  +  Tr{x,Xr)  =  i(0, 
dXr  _  dt 


(2.10) 

(2.11) 


2.2  Solutions  for  Constant  Velocity 

In  particular,  suppose  that  the  migration  velocity  is  a  constant,  c,  and  the  datum 
surface  is  the  a;-axis,  that  is,  Xs  =  (x3,0)  and  Xr  =  (xr,0).  Then 

T^{Xs,x)  =  pjc,  Tr{x,Xr)  =  pr/c,  (2.12) 

where 

ps  =  yJ{Xs- xY -\r  Pr  =  y/{Xr-xY  +  z^. 

For  this  case,  equation  (2.5)  and  (2.6)  become 


Ps+Pr-  Ct{(), 

(2.13) 

dps  dpr  dt 

^  “""dr 

(2.14) 

For  the  common-offset  case,  ^  =  y.  Then  differentiating  here  with  respect  to  y  and 
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using  equation  (2.7)  yields 


d/Jg  Xg  X  d  Py  Xy  X 

dy  Pi  ’  dy  py  ' 

Thus,  equations  (2.8)  and  (2.9)  are  simplified  to  algebraic  forms: 


+  /5r  =  ct{y,h), 


(2.15) 


Xs  —  X 


Ps 


+ 


Xy  —  X 


=  c 


dt 

dy' 


(2.16) 


Although  the  above  equations  are  still  not  easy  to  solve,  there  is  an  explicit  solution 
for  the  zero  offset.  In  fact,  h  =  0  implies  that  Xg  =  Xy  =  y  and  pg  =  p^  z=  sj{y  —  xY  +  z^. 
Therefore,  equations  (2.15)  and  (2.16)  are  further  simplified  to 


yj {y  -  xY =  |t(y), 
y  —  x  _  cdt 

^J{y  —  xY  +  ^ 

These  two  equations  yield  the  solution 


(2.17) 

(2.18) 


c‘  .  .  dt 


and 
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In  general,  equations  (2.8)  and  (2.9)  cannot  be  solved  explicitly.  However,  these 
two  equations  represent  a  general,  quantitative  relationship  between  migration  imaging 
and  velocity  distribution.  This  relationship  will  be  helpful  to  derive  formulas  for  velocity 
estimation. 


2.3  3-D  Case 

For  the  3-D  case,  the  derivations  and  formulas  are  completely  analogous  to  the  2-D 
case,  so  I  omit  a  detailed  derivation  and  just  list  the  main  results. 

In  the  3-D  case,  two  parameters  are  required  to  describe  the  datum  surface:  ^  = 
(^1)^2))  6 — iu-line  direction  and  ^2 — cross-line  direction.  Parallel  to  the  definitions  in 
the  2-D  case,  I  define 

X  =  (xi,a;2,2), 

Xy  —  (3rir(^),  :E2r(£))  ■2^r(^))- 
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Unless  specified  otherwise,  symbols  with  index  1  denote  the  variables  along  the  in-line 
direction  and  symbols  with  index  2  denote  the  variables  along  the  cross-line  direction. 
Snell’s  law  is  stated  as 


[Va;  r,(a;„  a;)  +  V®  Tr{x,  JCr)]  ’  ^  =  0’  *'  =  1,2. 

Now,  the  imaging  equations  become 

T^{Xs,x)  +Tr{x,Xr)  =  t{^), 
dTs{x„x)  dTr{x,Xr)  _  dt  ^ 

Let  h  =  (hi,  ^2)  denote  the  offset  vector  and  y  =  (yi,  2/2)  denote  the  midpoint  vector. 
In  the  common-ofiset  case,  the  imaging  equations  are 


(2.19) 

(2.20) 
(2.21) 


Ts{x^,x)+Tr{x,Xr)  =  t{y,h), 

£Zi  +  £i  =  ii  i  =  i  2 
%  ay,’ 


(2.22) 

(2.23) 


For  the  zero  offset  and  a  constant  velocity,  c,  there  is  an  explicit  solution  for  the  above 
equations: 


x  =  y--t{y)Vytiy), 


and 


2.4  Summary 

The  imaging  equations  provide  a  general  kinematic  representation  of  migration  pro¬ 
cesses.  Some  applications  of  these  equations  to  residual  migration  have  been  addressed  in 
the  literature  (Etgen,  1989;  Zhang,  1991).  In  this  thesis,  I  use  these  equations  to  formu¬ 
late  algorithms  for  migration  velocity  analysis.  From  the  derivation  of  these  equations,  it 
is  clear  and  meaningful  that  the  envelope  condition  (2.6)  and  Snell’s  law  (2.2)  are  equiv¬ 
alent.  Because  Snell’s  law  is  a  general  statement  in  reflection  seismology,  it  is  concluded 
that  the  imaging  equations  hold  for  converted  waves  in  isotropic  or  anisotropic  media, 
as  well.  That  means,  velocity  analysis  algorithms  derived  from  the  imaging  equations  in 
acoustic  media  can  be  generalized  to  converted  waves  and  anisotropic  media. 
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Chapter  3 


IMAGE  ANALYSIS  IN  A  CIG 


Velocity  analysis  must  be  done  on  multichannel  data.  After  prestack  migration,  the 
imaged  traces  are  sorted  into  common  image  gathers  (CIGs).  The  imaged  depths  in  CIGs 
provide  information  on  migration  velocity.  In  this  chapter,  I  am  going  to  quantify  the 
relationship  between  the  imaged  depths  and  the  velocity,  based  on  the  imaging  equations. 
For  simplicity,  suppose  that  the  migration  velocity  is  a  constant,  c,  and  the  datum  surface 
is  the  x-axis. 

By  the  definition  of  a  CIG,  to  study  images  in  a  CIG,  one  should  fix  the  horizontal 
coordinate,  x,  of  the  reflection  point.  Then,  both  the  imaged  depth  z  and  the  source- 
receiver  position  parameter  ^  can  be  considered  to  be  functions  of  the  velocity  c.  That 
is,  when  the  velocity  changes,  the  migration  image  at  a  particular  trace  location  changes 
in  the  vertical  direction.  Therefore,  z  is  a  function  of  c.  Similarly,  when  one  keeps  x 
unchanged,  the  image  at  x  will  come  from  difi'erent  source-receiver  pairs  for  dilferent 
choices  of  c.  This  fact  implies,  ^  is  a  function  of  c,  too.  Differentiating  equation  (2.13) 
with  respect  to  c. 


dps  ^1  ^ 

dz  dz  dc 


dps  dpr 

ae  ^  d^ 


dc 


dt  d^ 


(3.1) 


Now,  by  using  equation  (2.14)  to  eliminate  dps/d^  and  dprjd^,  and  using 


dps  _  Z  dpr  _  z 

dz  Ps  '  dz  Pr' 

I  obtain 

{Zjps  +  Zjpr)  —  =  t(0  =  {ps  +  Pr)/C. 

(3.2) 

(3.3) 

Solving  for  dz/dc,  then. 

dz  ^  P.P,  ^ 
dc  cz 

(3.4) 

Introduce  angles  9  and  0  as  in  Figure  3.1.  Then, 

z  z 

cos{9  —  </>)’  cos{6  +  4)  ’ 

(3.5) 

dz  z 

dc  ccos{6  —  (p)  cos{9  +  (j))' 

(3.6) 

Equation  (3.4)  shows  that  the  imaged  depth  increases  with  the  migration  velocity 
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Fig.  3.1.  Sketch  of  ray  paths,  (j)  is  the  dipping  angle  of  the  reflector,  and  6  is  half  the 
angle  between  the  source  and  receiver  ray  paths. 


for  fixed  x,  no  matter  what  the  true  medium  velocity  is.  This  conclusion  is  true  even  for  a 
general  migration  velocity.  Further  discussion  can  be  found  in  Appendix  A,  where  I  will 
show  this  result  for  downward  propagating  waves  and  also  show  that  the  sign  is  opposite 
for  upgoing  waves.  Moreover,  equation  (3.4)  allows  us  to  do  quantitative  analysis,  so 
one  may  find  how  the  degree  of  the  image  distortion  depends  on  the  magnitude  of  the 
velocity  error  and  on  the  refiection  geometry. 

Suppose  that  there  are  multichannel  data.  For  a  given  migration  velocity,  c,  common 
image  gathers  are  sorted  from  the  migrated  data.  Let  c*  denote  the  true  medium  velocity. 
At  c  =  c*,  the  imaged  depths  in  a  CIG  must  be  the  same,  i.e.,  z  is  independent  of  pairs 
of  {xa,Xr}  for  a  fixed  x.  With  no  loss  of  generality,  I  suppose  that  the  source  is  to  the 
left  of  the  receiver,  that  is,  0  >  0. 

3.1  Common-shot  Case 

In  the  common-shot  case,  different  imaged  traces  in  a  CIG  correspond  to  different 
shot  positions  at  a  fixed  trace  location.  Thus,  our  objective  here  is  to  describe  the 
sensitivity  of  the  velocity  analysis  technique  to  changes  in  shot  location  over  the  gather 
of  shots  that  image  a  particular  reflector  at  the  given  trace  location.  From  Figure  3.1,  I 
have 

Xg  —  X  = —ztan(0  —  (f>).  (3.7) 

As  in  the  discussion  above,  at  c  =  c*,  the  true  velocity,  the  imaged  depth  in  a  CIG  is 
independent  of  shot  position;  that  is. 


dz 

dxg 


C=C* 


=  0. 


(3.8) 
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If  one  moves  from  one  shot  position  to  another,  i.e.,  differentiating  with  respect  to  in 
equation  (3.7)  and  evaluating  the  result  at  c  =  c*  yields 


cos2(0  -  <f))  dxs  .  ’ 


dO  _  cos‘^{9  —  (f)) 

dx.  z 


(3.10) 


By  differentiating  equation  (3.6)  with  respect  to  Xs  and  using  equations  (3.8)  and  (3.10) 
I  find  that 


sin  26 


(3.11) 


dXadc\^_^,  CCOS^{6  +  4))  ’  V 

This  tells  us  that  the  depth  deviation  decreases  as  the  source  moves  right  {Axg  >  0) . 

Now  I  use  (3.11)  to  obtain  error  estimates  in  velocity  analysis.  From  Figure  3.1  I 
find  that 

sin  2^  _  cos{9  +  4>)  ^ 


X^  —  X. 


cos{9  +  (j))  =  — . 

Pr 


(3.12) 


The  first  equation  here  is  a  result  of  the  law  of  sines;  the  second  arises  from  the  trigonom¬ 
etry  of  the  right  triangle.  Furthermore,  Xj.  is  the  specular  receiver  position  corresponding 
to  the  shot  Xa  and,  hence  is  a  function  of  Xa  in  the  CIG.  Now  equation  (3.11)  becomes 


d'^Z  _  {Xa  -  Xr)  Pr 

dXadc  .  C’^Zpa 


(3.13) 


This  result  and  (3.8)  can  now  be  used  to  obtain  an  approximation  for  the  slope  of  the 
image  position  with  respect  to  changes  in  Xa  for  c  near  c*.  The  result  is 


dz  dz  d^z 

dXa  dXa  ,  dXa  dc 

a  a  C=C*  ^  C—C* 


(c  -  c*) 


{^a  Xj.)  Pri^C  C  ) 


C*  Z  Pa 


(3.14) 


Suppose  that  there  are  two  shots,  Xa^^Xa^,  and  <  Xa^.  Then  the  difference  in  imaged 
depths  between  these  two  shots  in  a  CIG  can  be  approximated  by 

z{Xa,)  -  z{Xa,)  «  (x,,  -  X, J  «  {Xa^  -  xj  (x,,  -  (3.I5) 

OXg  C  Z  pg 

In  this  equation,  Xa^  =  (x^^  -f  Xa^)l2  and  Xrp  is  the  corresponding  receiver  position 
producing  the  specular  image  of  interest  on  the  trace  at  x.  Thus,  solving  for  the  relative 
velocity  error,  I  find  that 


Pa  _  AZZ  Pa 

C  (^^0  ^ro)  (^52  ^51 )  Pv  (^30  ^tq)  pr 


(3.16) 


Here,  Az  =  z(Xaj)  —  z(Xjj)  and  Ax^  =  x^^  —  x^, .  This  relationship  is  valid  when  c 


15 


approximates  c*  and  the  two  shots  are  close  for  each  other.  The  quotient  is  greater 
than  one  for  the  negative  dip  angle,  and  less  than  one  for  the  positive  dip  angle. 

In  equation  (3.16)  one  see  the  factors  that  govern  the  accuracy  of  velocity  analysis, 
(c  —  c*)/c*,  under  the  assumption  that  the  medium  velocity  is  constant.  The  accuracy 
of  velocity  analysis  for  a  CIG  of  common  shot  images  is  best  for  large  source-to-receiver 
offsets,  well-separated  shot  points,  and  shallower  targets,  which  is  well-known  but  not 
previously  quantified.  Interestingly,  it  is  also  better  for  reflectors  with  positive  dip  (i.e., 
receivers  located  in  the  down-dip  direction  relative  to  the  shot  point). 


3.2  Common-offset  Case 


In  the  common-offset  case,  different  imaged  traces  in  a  CIG  correspond  to  different 
offsets.  Thus,  I  seek  a  description  of  the  dependence  of  the  image  depth  2:  on  the  half¬ 
offset  h  at  a  fixed  trace  location,  x.  From  Figure  3.1, 


2h  =  Xr  —  Xa  =  2(.tan(0  —  <f>)  +  tan(^  -f-  0)). 


Similar  to  the  derivation  of  (3.11),  I  find  that 


d^z  _  2  sin  26 

dhdc  _  ,  c  {cos^{6  +  (/))+  cos^{9  —  d>)) 


(3.17) 


(3.18) 


This  result  tells  us  that  the  depth  deviation  increases  as  the  offset  increases  (A/i  >  0). 
By  using  the  relations  in  equation  (3.12),  equation  (3.18)  becomes 


d'^z 

dhdc 


2  h  2  PaPf. 


(3.19) 


dz  2 /i  2  Pspr  ( 


(3.20) 


Suppose  that  there  are  two  offsets  hi,  7*2  and  hi  <  /i2-  Then  the  difference  in  imaged 
depths  for  the  two  offsets  in  a  CIG  is  given  by  the  approximation 


z(h2)  -  z{hi)  w  — (h2  -  hi)  «  2  ho  (h2  -  hi)  ^  , 

Oh  c*z  p2  +  p2 

where  ho  =  (h2  -I-  hi)/2.  Thus, 

(c  -  c*)  {z{h2)-z{hx))zpl-\-pl  ^  Azz  p^Pr 

C*  2  /iQ  (^2  ~  ^1)  ^PsPr  2  Hq  ^PsPr 


(3.21) 


(3.22) 


where  Az  =  z(h2)  —  -^(hi)  and  Ah  =  h2  —  hi.  This  relationship  is  valid  when  c  approxi¬ 
mates  c*  and  the  two  offsets  are  close  for  each  other.  The  quotient  {pi  -f  pi) l2paPr  equals 
one  for  a  horizontal  reflector  and  is  greater  than  one  for  any  given  dip  angle. 
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The  relationship  (3.22)  shows  us  that  the  accuracy  of  velocity  analysis  for  a  CIG  of 
common  offset  images  is  best  when  the  two  offsets  are  large  and  greatly  different  from  one 
another,  and  shallower  targets,  which,  again,  is  well-known  but  not  previously  quantified. 
Interestingly,  this  accuracy  deteriorates  with  increasing  dip. 

Note.  From  the  result  (3.22),  one  can  conclude  that  any  error  of  velocity  and  the 
difference  of  the  offsets  results  in  a  nonzero  deviation  Az.  More  precisely, 

c=c* 

Therefore,  I  define  the  quantity  d^zfdhdc  as  the  sensitivity  to  the  velocity  error  for 
common  offset  images.  Similarity,  the  quantity  d^zfdxsdc  can  be  defined  as  the  sensitivity 
to  the  velocity  error  for  common  shot  images. 

3.3  Multiple-layer  Case 

If  the  medium  is  made  up  of  more  than  one  layer,  the  expression  for  the  error  estimate 
should  be  modified.  Let  us  consider  only  a  simple  model  that  consists  of  two  horizontal 
layers,  in  the  common-offset  situation  as  in  Figure  3.2.  Differentiating  equation  (2.1) 


Fig.  3.2.  Two-layer  model.  ci  and  C2  are  layer  velocities,  di  and  d2  are  layer 
thicknesses,  and  0i  and  62  are  angles  of  raypath  from  the  vertical. 


with  respect  to  ca,  and  using  equation  (2.2)  yield 


dTs{Xs,x)  dTr{x,Xr) 


dz 


dz 


dz 

dC2 


dTs{Xs,x)  dTr{x,Xr) 


dC2 


dco 


By  using  the  notations  shown  in  Figure  3.2,  the  traveltimes  are  represented  by 

di  _  d2 


T^{x^,x)  =  Tr{x,Xr)  = 
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Cl  cos  01  C2COS02 


(3.23) 


(3.24) 


Although  6 1  and  62  are  functions  of  C2,  two  variations  from  61  and  02  will  be  balanced 
to  each  other  from  Fermat’s  principle.  Therefore,  differentiating  equation  (3.24)  with 
respect  to  C2  will  yield 


dts  _  dTr  _  d2 

(3.25) 

dc2  dc2  cos  $2 

Substituting  this  result  and 

dTg  dTr  cos  02 

dz  dz  C2  ' 

(3.26) 

into  equation  (3.23),  I  obtain 

dz  d2 

dC2  C2  COs2  02  ' 

(3.27) 

Now  I  want  to  derive  a  formula  for  the  sensitivity  d‘^zldc2dh.  Firstly,  I  evaluate  C2  at  the 
true  velocity  of  the  second  layer,  so  that  the  imaged  depth  z  and  <^2  will  be  independent 
of  offset.  Introduce  the  horizontal-slowness  parameter 


Then 


sin^i  _  sin02 

Cl  C2 


cos  01  =  \/l  —  Cip2, 


cos 


02=\/r 


/•2„2 

C2P  7 


(3.28) 

(3.29) 


h  = 


dicip 


+ 


d2C2P 


dp 


ijl  -  cfp2  yj\  -(^2 
diC\  d2C2 


+ 


(1  -  cfp2)3/2  ^  (1  _  f^2)3/2- 

From  equations  (3.27)  and  (3.29),  I  have 

d^z 
dhdc2 

By  using  equation  (3.31),  the  above  equation  becomes 


fk  A  j 

f  ’  1 

1  _  j 

(  ^  ] 

j  dp  2d2C2P  dp 

C2  dh  ' 

i,COs2  02  J 

C2  dp  ' 

'  dh  (1  —  C2p2)2  dh 

d‘^z 


2d2C2P 


dhdC2  (1  -  cip2)2 


diCi 


d2C2 


T-1 


(1  -  C?p2)3/2  ^  (1  _  ^2)3/2j 


When  C2  «  Cl,  equation  (3.33)  is  simplified  to 


(3.30) 

(3.31) 


(3.32) 


(3.33) 


d'^z  _  d2  2  tan  0i 

dhdc2  di+d2  Cl 


(3.34) 


Compared  to  equation  (3.18)  (let  0  =  0),  equation  (3.34)  shows  us  that  for  multiple 
layers,  the  sensitivity  to  the  velocity  error  in  a  thin  layer  is  reduced  by  the  ratio  of  the 
layer  thickness  to  the  reflector  depth. 
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To  go  a  further,  I  introduce  the  ratio  of  the  layer  velocities,  77  =  C2/C1,  then  equa¬ 
tion  (3.33)  becomes 


d'^z 


2d2r]p 


-I- 


d^T) 


-1 


dhdc2  (1  —  [(1  —  (1  — 

When  T]  is  small,  the  above  equation  can  be  approximated  by 


(3.35) 


d^z  o  do 

—  «2wcose,-.  (3.36) 

Equation  (3.36)  quantifies  how  the  sensitivity  to  the  velocity  error  deteriorates  for  a  low- 
velocity  zone  of  the  target  layer. 

Moreover,  for  n  horizontal  layers,  equation  (3.33)  can  be  generalized  to 


dhdc„ 


2djiCjiP 

(1  -  C2p2)2 


dfCf 

(1  —  cfp2^3/2  ’ 


(3.37) 


where  c,-  is  the  7-th  layer  velocity  and  d,-  is  the  thickness  of  the  i-th  layer. 

For  a  more  complex  velocity  model,  the  formula  for  the  sensitivity  of  velocity  analysis 
could  not  be  so  explicit,  but  it  is  still  computable.  The  detailed  discussion  will  be  seen 
in  Section  5.2. 
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Chapter  4 


HYPERBOLIC  RESIDUAL  MOVEOUT 


When  an  incorrect  velocity  is  used  to  migrate  multichannel  data,  the  imaged  depths 
in  a  CIG  will  differ  from  each  other.  In  this  situation,  residual  moveout  (RMO)  is 
observed  in  migrated  data.  Like  normal  moveout,  residual  moveout  contains  information 
from  which  one  can  estimate  the  medium  velocity.  Some  geophysicists  have  applied  this 
technique  to  migration  velocity  analysis.  They  assume  a  hyperbola  for  residual  moveout 
and  estimate  an  average  velocity,  then  convert  this  average  velocity  into  the  interval 
velocity.  In  this  chapter,  I  develop  a  theoretical  basis  for  this  kind  of  application.  Under 
the  assumption  of  small  offsets,  I  derive  general  representations  for  residual  moveout  by 
using  the  imaging  equations  introduced  in  Chapter  2.  Based  on  these  formulas,  I  define 
RMO  velocity  and  show  that  residual  moveout  is  dominated  by  the  difference  between  the 
RMO  velocity  and  migration  velocity.  For  a  laterally  invariant  velocity,  I  prove  that  the 
RMO  velocity  is  a  good  approximation  to  the  rms  velocity,  so  one  can  directly  estimate 
velocity  by  residual-moveout  correction.  In  this  way,  velocity  analysis  requires  only  a 
single  prestack  migration,  so  one  can  avoid  relatively  costly  remigrations.  When  the 
velocity  has  a  weak  lateral  anomaly  and  the  reflector  is  horizontal,  I  develop  a  formula 
to  calculate  the  interval  velocity  from  the  RMO  velocity.  For  a  complex  structure,  I  also 
show  the  limitation  of  this  approach  for  velocity  estimation. 

4.1  Residual  Moveout  Representations 

Consider  the  2-D,  common  offset  case.  All  notations  have  the  same  meaning  as 
those  in  Chapter  2.  Suppose  that  prestack  migration  is  implemented  on  a  multi-offset 
data  set  by  using  an  initial  velocity.  The  migrated  data  is  sorted  into  CIGs.  At  each 
CIG,  the  imaged  depth  is  a  function  of  offset,  i.e.,  z  =  z{h)  (see  Figure  4.1).  When  the 
migration  velocity  is  correct,  z{h)  =  z(0)  for  all  offsets.  Otherwise,  for  a  wrong  velocity, 
one  should  expect  that  z{h)  ^  -^(0).  It  is  desirable  to  know  the  relationship  between 
residual  moveout,  z{h)  —  ^(0),  and  the  error  in  velocity.  Below,  I  derive  this  relationship 
for  the  simplest  case  and  then  generalize  that  result. 

4.1.1  Constant  velocity  and  horizontal  reflector 

Suppose  that  the  medium  velocity  is  a  constant  v,  migration  velocity  is  a  constant 
c,  and  the  reflector  is  horizontal.  Then  the  recorded  traveltime  is  given  by 

t‘^{y,h)  =t‘^{y,0) +  ^,  (4.1) 
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z 

Fig.  4.1.  a  Sketch  for  Common  Image  Gather. 

while  the  traveltime  to  the  migration  image  is  given  by 

t{y,  /.)  =  2^  =  ^  =  H  VFT7?.  (4.2) 

Solving  here  for  the  image  depth  2  as  a  function  of  offset  yields 

h)-h‘^  =  0)  +  =  ^2(0)  +  /i2  ^4  3^ 

That  is,  the  moveout  is  a  hyperbola  for  c  >  t;  or  an  ellipse  for  c  <  v: 

z\h)  =  z2(0)  +  -  l)  (4.4) 

Equation  (4.4)  shows  that  the  residual  moveout  is  the  exact  hyperbola  or  ellipse  for 
constant  velocities  v  and  c,  and  horizontal  reflector.  This  result  is  parallel  to  that  of 
moveout  in  unmigrated  data. 

4.1.2  General  case 

For  a  general  velocity  or  an  arbitrary  reflector,  one  should  not  expect  a  simple 
expression  such  as  (4.4).  However,  the  more  general  result  can  be  approximately  derived 
by  considering  an  asymptotic  expansion  for  small  offset.  Thus,  I  seek  a  Taylor  series  for 
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z  as  a  function  of  h. 

First,  because  2  is  a  symmetric  function  of  h,  1  obtain 


dz_  S 

^  /.=o  ^  ^ 

This  result  implies  the  Taylor  series  expansion 


0,  etc. 


z^{h)  —  2^(0)  +  2~^j^  0{h‘^).  (4.6) 

/i=0 

Now  let  us  try  to  estimate  the  second  derivative  term. 

For  fixed  x,  the  midpoint  y  and  the  imaged  depth  2  are  functions  of  offset  h.  Dif¬ 
ferentiating  equation  (2.8)  with  respect  to  h,  yields 


dr^  ,  dTr]  dy  ,  ,  dr, 


dy  dy\  dh  ^  dh 


-1  \dr, 
'\  [dz 


dz  dt  dy  dt 


dz  \  dh  dy  dh  dh' 


Using  equation  (2.9)  to  cancel  the  first  terms  on  left  and  right,  I  obtain 


dxs  dtr  dz  _  dt  dxa  dr/ 

dz  dz  dh  dh  dh  dh 


Notice  that  y  is  symmetric  in  h,  so 


dy 

=  0. 


Using  equations  (4.5)  and  (4.9),  I  find  the  derivatives  of  both  sides  in  equation  (4.8), 
with  respect  to  h  at  /i  =  0,  are  given  by 


d  dt  f  dr^ 
dh  dh  i  dh 


Srr'Nl  dH  r^^'Ts  . 


d  /  dxg  dxr  dz 

dh  I  ^2  ^2  j  dh 


dr,  dTr\  d?z 


dz  dz  I  dh'^ 


(4.11) 


Differentiating  equation  (4.8)  with  respect  to  h,  and  using  (4.10)  and  (4.11),  I  set  up  the 
following  equation  for  (fz/dh?  |/i=o: 


dr^  ^ 

dz  dz  }  dh^ 


d'^T,  d\ 

dh‘^  a/i2 


(4.12) 


Equation  (4.12)  holds  for  any  medium  velocity  function,  any  migration  velocity,  and  an 
arbitrary  reflector.  Now  I  will  simplify  this  equation  for  the  case  of  constant  migration 
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velocity. 

Suppose  the  migration  velocity  is  a  constant,  c;  then 

Ts  =  \Jz^  ->r{x-y  +  h)2/c,  Tr  =  \J  {x  -  y  -  hy/c, 

and, 

Bts  _  dxr  _  z 
h=o  ~  h=o  ~ 

d'^Ts  _  d'^Tr  _  1  (a:  —  yY 

dh?  ^^0  “  dh?  “  cp  cp^ 

where 

p  =  sjz‘^->r{x-  y)2. 

Equations  (4.15)  and  (4.16)  yield 

dTa  dXr  _  2Z 
dz  ^  ”  cp’ 

d'^Ts  d'^Tr  _  2  2(a:  —  yY 

dh"^  dh"^  cp  cp^ 

Using  these  results,  equation  (4.12)  is  simplified  by 

2z  (fz  _  dH  2{x  —  yY  2 

.(^pdh‘^\h=o~  h=o'^  cp3  cp' 

Furthermore, 

\  d?z]  _ 

~  dh^  ’ 

L  J/i=o  h=0 

c 

p  =  -t  \h=0, 

{y  —  x)  c  dt 

P  ~  2  dy 

Thus,  I  obtain  the  result 


(fz‘^  _c‘^\dH 

dh2  “  2  ^dh?  ^  V^j 
'»=o  L  >  U=o 

(PtI  -o\  ^  f— VI  -- 

dh?  .  „  ^  dh?  ”*"15?/]  c2 

'^=0  L  ^  J/i=o 


(4.13) 

(4.14) 

(4.15) 

(4.16) 

(4.17) 

(4.18) 

(4.19) 

(4.20) 

(4.21) 

(4.22) 

(4.23) 

(4.24) 
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where  is  the  migration  time  defined  by 


Tm  = 


C  ’ 


(4.25) 


From  equation  (4.23),  it  is  concluded  that  after  migration,  the  the  main  part  of  the 
residual  moveout  is  determined  by  the  migration  velocity  c,  curvature  of  the  unmigrated 
data  tdH/dh?,  and  slope  of  the  unmigrated  data  dt/dy.  This  statement  is  true  for  a 
constant  migration  velocity,  any  medium  velocity  distribution,  and  an  arbitrary  reflector. 
Based  on  equation  (4.23),  I  define  a  residual-moveout  velocity  (RMO  velocity)  by 


1  _  1  r  dH 

[vrmoY  4  dh?’'^\dy) 


By  using  the  RMO  velocity,  the  residual  moveout  is  written  as 


(4.26) 


’i('!)  =  ’i(0)+(-^-i)  4/i’  +  0(A'‘).  (4.27) 

\  ^rmo  ^  / 

Equation  (4.27)  shows  the  RMO  velocity  can  be  directly  estimated  from  the  residual 
moveout.  Compared  to  NMO  velocity,  a  new  term,  dt/dy,  is  involved  in  the  definition 
of  RMO  velocity,  and  these  two  velocities  will  be  equal  if  and  only  if  dt/dy  =  0.  In 
the  following  examples,  it  is  seen  that  this  term  removes  the  dip  elFect  in  RMO  velocity 
analysis  and  makes  the  RMO  velocity  approximately  equal  to  rms  velocity  under  the 
condition  of  lateral  velocity  homogeneity.  In  addition,  the  RMO  velocity  analysis  based 
on  equation  (4.27)  is  actually  composed  of  an  inverse  NMO  with  velocity  c  and  an 
NMO  velocity  analysis.  Therefore,  all  programs  for  NMO  velocity  analysis  can  be  used 
to  do  RMO  velocity  analysis.  This  technique  was  applied  by  Deregowsky  (1990).  In 
his  paper,  Deregowsky  used  an  iterative  algorithm  in  RMO  velocity  analysis  to  handle 
lateral  velocity  variation.  However,  this  iteration  may  fail  as  I  will  show  later. 


4.2  Velocity  Estimates  for  Constant  Migration  Velocity 

For  arbitrary  migration  velocity,  the  residual  moveout  representation  (4.12)  will  be 
too  complex  to  yield  a  quantitative  formula  for  velocity  estimates.  In  contrast,  when 
migration  velocity  is  a  constant,  the  residual  moveout  representation  is  a  more  explicit 
function  of  the  velocity  error,  which  allows  velocity  estimates  to  be  done  without  iteration. 
Furthermore,  prestack  migration  with  a  constant  velocity  is  most  efficient. 

In  addition  to  the  normal  moveout  formula  (4.4),  I  will  compute  residual  moveout 
by  the  formula  (4.23)  and  (4.24)  for  several  special  cases.  These  results  are  similar  to 
normal  moveout,  except  that  velocity  estimates  by  the  former  are  insensitive  to  reflector 
dips. 
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4.2.1  Constant  medium  velocity 
Dipping  reflector 

Suppose  that  the  medium  velocity  is  a  constant  u,  and  the  reflector  dip  is  6.  In  this 


case. 

t^{y,h)  =  Ay^s\v?9lv^  Ah?  co^  0  j . 

(4.28) 

Therefore, 

<0.2  =4cos^9/t)^, 

(4,29) 

dt 

—  ^2sin6/v, 

9y  h=o 

(4.30) 

From  equation  (4.23),  I  obtain 

(Pz^  /  9  /  9 

2  =  h  - 1). 

/j=0 

(4.31) 

z‘^{h)  =  ^2(0)  +  (c2/t;2  -  1)  h2  +  0{h?). 

(4.32) 

or 

Tl{h)  =  rl{Q)  +  {llv‘^-llc^)Ah'^  +  0{h% 

(4.33) 

Equation  (4.33)  shows  that  for  a  constant  medium  velocity  and  a  dipping  reflector, 
the  residual  moveout  is  independent  of  the  reflector  dip  and  the  RMO  velocity  equals  the 
medium  velocity. 


Diffraction  from  a  scattering  point 

Suppose  that  the  medium  velocity  is  a  constant  v,  and  a  scatterer  is  located  at  the 
point  (x*,z*).  In  this  case, 


z^(h)  =  z^(0)  +  -  1)  +  0(h'), 


(4.34) 

(4.35) 

(4.36) 

(4.37) 

(4.38) 
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or 

r^ih)  =  r2  (0)  +  (l/t;2  -  l/c^)  4/1^  +  (4.39) 

Again,  equation  (4.39)  shows  that  the  residual  moveout  is  independent  of  the  lateral 
offset  from  the  point  scatterer  and  the  RMO  velocity  equals  the  medium  velocity,  when 
the  medium  velocity  is  a  constant. 


4.2.2  Laterally  invariant  medium  velocity 
Horizontal  reflector 

Suppose  that  the  medium  velocity  is  a  laterally  invariant  function  v{z),  and  the 
reflector  is  horizontal  with  depth  z*.  In  this  case, 

t^{y,h)  =  t^{y,Q)  +  4h^ /v%^^{z*)  +  Oih"^),  (4.40) 


Where  Vrms  is  the  rms  velocity  of  v.  Differentiation  of  the  above  equation  yields 


h=Q 


Therefore,  from  equation  (4.23), 


=  0. 


(4.41) 

(4.42) 


(Pz^ 


dh^ 


h=0 


2P 


(4.43) 


and 


(4.44) 

^rmo('2^(0))  —  )• 

(4.45) 

Notice  that  2r(0)  ^  z*.  That  is,  in  depth  migration,  the  imaged  depth  is  inconsistent 
with  the  desired  point  at  which  the  rms  velocity  is  determined  from  the  residual  moveout. 
However,  if  let 


r*  =  2 


ds 

t;(s)’ 


(4.46) 


then 


ri{h)=Ti{0)  + 


4*2  +  c)(h'‘). 


(4.47) 


and  r,„(0)  =  r*.  Therefore,  time  migration  can  give  us  the  correct  location  at  which  the 
rms  velocity  is  determined  by  the  residual  moveout. 

Equation  (4.47)  shows  that  when  the  true  velocity  is  laterally  invariant,  the  RMO 
velocity  equals  rms  velocity  at  the  migration  time. 
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Dipping  reflector 

Suppose  that  the  medium  velocity  is  a  laterally  invariant  function  u(r),  and  the 
reflector  angle  is  In  this  case, 


dH 

dh^ 


dt 


dy 


=  2p, 


(4.48) 

(4.49) 

(4.50) 


/i=0 


where  r*  is  the  vertical  time  at  the  reflection  point  and  p  is  the  horizontal  slowness,  i.e., 

sin^ 


P 


Therefore, 


dh?  I  5?/ 


=  4 


h=Q 


v{t*) ' 

/o-*(l-pV)-3/2da 

Jf  u2(l  -pH‘^)-3W 


and 


Jq  v^{a){l  —  p'^v‘^{a)) 


/o’’*(l-pV((T))-3/2d(7  ■ 

Notice  that  Vrmo  does  not  equal  the  rms  velocity  that  is  deflned  by 


=  v^{<y)da. 


In  fact,  from  the  Taylor  series  expansion,  derived  in  Appendix  A, 


where 


vLo  =  vl  +  ^-{v\-vt)p-^  +  0{p\ 


v\  =  —  /  v'^{s)ds. 

T*  Jo 


(4.51) 


(4.52) 


(4.53) 


(4.54) 


(4.55) 


(4.56) 


V4  is  always  greater  than  V2  and  they  are  equal  only  if  v{z)  is  a  constant.  Moveover,  Vrmo 
is  close  to  V2{t*)  for  small  p  or  a  small  gradient  of  v{z).  When  Vrmo  is  estimated  from  the 
residual  moveout  formula  (4.27),  Vrmo  is  measured  at  the  migration  time,  rm(0)  instead 
of  r*.  That  means,  it  is  desirable  if 

^rmo  ~  ^2(An(0))-  (4.57) 

Fortunately,  Vrmo  is  a  better  approximation  of  U2('bn(0))  than  V2{t*).  In  fact,  if  is  a 
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linear  function  of  migration  time,  one  can  prove  that  v^mo  is  exactly  equal  to  U2('^m(0)) 
for  the  choice  c  =  U2(0),  no  matter  what  dip  the  reflector. 

It  is  concluded  that  for  choosing  c  =  ^^2(0)1  the  RMO  velocity  is  a  good  approxima¬ 
tion  of  the  rms  velocity  at  the  migration  time  as  long  as  the  velocity  varies  mildly. 


4.2.3  Higher  residual  terms 

The  residual  moveout  is  hyperbola-like  only  for  small  offset.  In  fact,  offsets  should 
not  be  too  small  so  that  resolution  in  velocity  analysis  is  not  too  low  (See  Chapter  3). 
Therefore,  one  should  consider  the  error  due  to  truncating  the  higher  residual  terms. 
This  work  is  partially  implemented  in  Mathematica. 

Using  the  rule  of  differentiation  for  a  compound  function  in  equation  (4.12)  and 
setting  c  constant,  the  fourth  order  derivative  of  z  with  respect  to  h  satisfles 


'd{p^+pr)  d^z' 

'  d^t  (fiy  d^t 

^{Ps  +  Pt) 

dz  dh^ 

—  c 

/i=0 

dh^  ^  dh^  dh^dy 

/i=0 

dh^ 

'(Pz  d^ip^+Pr)' 

'(Pyd^ip^  +  prY 

dh^  dh^dz 

—  0 

dh?'  dh'^dy 

(4.58) 


For  a  constant  velocity  v  and  a  dipping  reflector  with  angle  8,  equation  (4.58)  is  simplifled 
to 


dh^ 


h=0 


24{l/v^-l/(^)  sm^29 

{z*Y 


(4.59) 


where  r„,  is  the  migration  time  and  z*  is  the  depth  of  the  reflection  point.  Therefore, 
from  (4.59)  and  (4.33),  a  more  accurate  expression  for  the  residual  moveout  is 


(l/t;2 


—  1  /c^)  sin^  29 


h^  +  0{h^). 


(4.60) 


For  a  constant  velocity  v  and  a  scattering  point  at  (a:*,  2*),  again 


dh^ 


h=Q 


_  24{\ I —  1 1 (?■)  s\v?  29 

-  (?p  ’ 


except  that  now 


9  =  arctan 


y-x 


In  comparison,  for  the  unmigrated  data,  the  higher  residual  term  is 


dh^ 


24sin2  29cos^9 


h=0 


v^(z*y 


(4.61) 


(4.62) 


Equations  (4.59)  and  (4.61)  show  that  when  the  true  velocity  is  a  constant,  small  higher- 
residual  terms  of  moveout  are  obtained  for  a  close  background  velocity,  small  offset,  small 
ratio  of  offset  to  the  imaged  depth,  and  the  dipping  angle  that  is  near  0  or  90  degrees. 
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For  a  laterally  invariant  velocity  and  a  horizontal  reflector  at  the  zero-ofiset  time  to, 


^  24{vj  -  v\) 

h=Q  ^0^2 


(4.63) 


where  V2  and  V4  are  defined  in  equations  (4.53)  and  (4.55). 

Equation  (4.63)  shows  that  when  the  medium  velocity  is  laterally  invariant,  a  small 
higher-residual  term  is  obtained  for  small  gradient  of  the  velocity. 


4.2.4  Lateral  velocity  anomaly 

RMO  velocity  can  be  estimated  from  residual  moveout.  If  the  medium  velocity  is 
laterally  invariant,  the  RMO  velocity  approximates  the  rms  velocity.  However,  when  the 
velocity  has  a  lateral  anomaly,  this  conclusion  fails  so  that  one  cannot  calculate  directly 
the  interval  velocity  from  the  RMO  velocity  by  Dix  equation  or  other  algorithms.  Now  I 
will  develop  an  equation  to  solve  for  the  laterally  variant  interval  velocity.  This  equation 
holds  for  a  small  lateral  velocity  variation  and  a  horizontal  reflector.  For  horizontal 
reflectors,  the  RMO  velocity  is  the  same  as  the  asymptotic  stacking  velocity  in  normal 
moveout. 

Suppose  the  true  slowness  w{x,  z)  can  be  written  as 

w(x,z}  =  w(z)(l+a(x,z)),  (4.64) 


where  iv(z)  is  a  reference  slowness  and  a(x,z)  is  a  small  perturbation.  The  slowness 
w{x,z)  and  horizontal  reflectors  will  generate  a  two-way  traveltime  map,  t{y,h).  After 
prestack  migration  with  an  initial  migration  velocity,  a  RMO  velocity  (therefore,  slow¬ 
ness)  can  be  calculated  from  residual  moveouts.  Let  Wa{x,z)  denote  this  RMO  slowness, 
Ws{z)  denote  the  average  of  w^{x,z)  over  the  a;-direction;  then  an  equation  for  a{x,z)  is 
obtained  in  Appendix  B: 


w^{z) 


- 1 


=  '/ 

tn  Jo 


da 

v{ay 


where  Vg  =  l/wg,  to  is  the  two-way  zero-offset  time,  and  v  be  the  solution  of 


(4.65) 


(4.66) 


When  a  and  w  are  independent  of  depth,  the  result  in  equation  (4.65)  is  the  same  as 
that  of  Lynn  and  Claerbout  (1982).  Equation  (4.65)  shows  that  the  second  derivative  of 
a  determines  the  anomaly  of  the  RMO  velocity;  the  anomaly  of  the  RMO  velocity  at  a 
depth  results  from  the  anomaly  of  the  interval  velocity  above  this  depth.  Furthermore, 
vds  increases  as  a  decreases,  so  the  anomaly  of  the  interval  velocity  near  the  surface 
has  the  largest  effect. 


29 


Applying  the  Fourier  transform,  with  respect  to  x,  to  equation  (4.65),  I  obtain 


2 

g{k^,z)  =  jr  a{kx,(7) 
60  *'0 


+  1  + 


V^{<7) 


da 

v{a)  ’ 


(4.67) 


where  kx  is  the  wavenumber  and  g  is  the  Fourier  transform  of  —  1.  Equation  (4.67) 

is  a  First-kind  Volterra  integral  equation  that  is  ill-posed.  Therefore,  the  recursive 
alogrithm  for  equation  (4.67)  is  unstable.  To  obtain  a  stable  solution,  one  may  apply  the 
damped  least-squares  method  to  this  equation. 


4.3  Some  Extensions 

The  previous  discussion  in  this  chapter  focused  on  the  common  offset  case  in  2-D. 
Now  I  will  discuss  some  extensions  to  common  shot  and  3-D  situations. 


4.3.1  Common  shot  case 

Although  common  offset  migration  and  common  shot  migration  are  stated  in  differ¬ 
ent  forms,  they  should  represent  the  same  physical  process  and  produce  the  same  images, 
if  the  same  multichannel  data  are  used.  Particularly,  for  a  constant  migration,  the  resid¬ 
ual  moveout  in  common  shot  should  have  the  same  representation  as  equation  (4.27) 

^m(^)  =  T-m(O)  +  -  4)  4^2  -h  0(/i^).  (4.68) 

The  problem  is:  for  common  shot,  migrated  data  in  a  CIG  are  sorted  according  to 
shot  positions  instead  of  the  shot-receiver  offsets.  That  means,  one  does  not  directly  know 
what  offset  values  the  migrated  traces  in  a  CIG  have,  unless  the  reflector  is  horizontal. 
In  order  to  overcome  this  difSculty,  the  offsets  must  be  calculated  for  the  traces  in  a  CIG. 
Lee  and  Zhang  (1991)  proposed  a  formula  to  implement  this  work.  However,  Lee  and 
Zhang’s  formula  requires  information  of  reflector  dip  and  holds  for  small-dip  reflector  and 
constant  medium  velocity.  Therefore,  the  application  of  this  formula  is  limited. 

Another  way  to  calculate  the  offset  is  by  using  the  Kirchhoff  integral.  In  the  Kirchhoff 
summation,  one  calculate  two  inversion  outputs  which  have  two  different  amplitudes.  One 
is  the  original  amplitude;  the  other  is  the  original  one  multiplied  by  offset  h.  Thus,  the 
ratio  of  the  amplitudes  of  the  two  outputs  will  give  the  offset  value  at  the  specular  source- 
receiver  position.  In  this  way,  a  mapping  between  shot  position  and  offset  is  obtained. 
The  same  technique  was  used  to  determine  the  angle  of  reflection  in  Kirchhoff  inversion. 
[See  Bleistein  et  al.  (1987).] 
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4.3.2  3-D  czise 


Consider  the  3-D  common  offset  case.  The  imaging  equations  are  stated  in  equa¬ 
tions  (2.22)  and  (2.23).  The  Taylor  series  of  the  imaged  depth  is  expanded  into 


z\h)  =  -K  ^ 


2  dh\ 


/i=0 


h\  + 


d^z^ 


dhidh^ 


^  ^  15V 

hih2  -f- 

/i=0  2  dhl 


/i=0 


hl  +  0{\h\^),  (4.69) 


where  the  odd-order  derivatives  are  eliminated  because  of  the  symmetric  condition 

z{h)  =  z{-h).  (4.70) 

Similar  to  the  derivation  of  equation  (4.12),  each  of  the  second  derivatives  in  equa¬ 
tion  (4.69)  will  individually  satisfy 


'dr.  dTr\  d‘^z 


+ 


dz  dz  J  dhidhjj 


dH 


/i=0  dhidhj 


/i=0 


^  d^Tr 


dhidhj  dhidhj] 


/i=0 


(4.71) 


where  the  indices,  i  and  j,  take  the  integer  value  1  or  2.  Equation  (4.71)  is  valid  for  any 
velocity  distribution  and  an  arbitrary  reflector. 

Suppose  that  the  migration  velocity  is  a  constant,  c;  then 


Ta  =  \fz‘^  +  \x  —  y  +  /ip/c,  Tr  =  /  c. 


So, 


and 


dn  ^ 

dz  dz 


-  ^ 
/i=0  cp’ 


5^r, 


-I- 


di^Tr 


dhidhj  dhidhj] 


/i=0 


_  2  ^  2(a:,'  yi){xj  yj) 


where  6ij  is  the  Kronecker  delta  (1  ii  i  =  j  and  0  if  zV  i),  and 


p  =  \Jz‘^  +  I®  -yp. 


Let  to(2/)  =  <(y,0);  then 


and 


P  =  2*0, 


[Vi  -  ^i)  _  cdto  .  .  „ 

o  A  ’  ^  1 7  2. 

p  2  dyi 


By  using  relations  (4.76)  and  (4.77),  equation  (4.74)  is  rewritten  as 


dhidhj  dhidhj 


4  \  dtQ  dto 

''On 


h=0  cHo  tudyidyf 


(4.72) 

(4.73) 

(4.74) 

(4.75) 

(4.76) 

(4.77) 

(4.78) 
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Furthermore,  suppose  that  the  medium  velocity  is  a  constant,  v.  If  c  =  v,  all  the 
derivatives  in  the  left  of  equation  (4.71)  will  be  zero  because  the  imaged  depth  z  is 
independent  of  offset  for  a  correct  migration  velocity.  That  is. 


dH 

dhidhj 


h=0 


d^T, 


+ 


d^Tr 


dhidhj  dhidhj 


h=0,c=v 


Using  equation  (4.78),  and  setting  c  =  u,  I  rewrite  equation  (4.79)  as 


(4.79) 


dH 


dhidhj 


4  1  ^0  dh 

~Vii 


h=0  ^^^0  todyidyj' 


(4.80) 


Since  t  is  just  the  recorded  traveltime,  the  above  expression  must  be  independent  of  the 
migration  velocity  c.  Substracting  equation  (4.78)  from  equation  (4.80),  yields 


dH 


dhidhj 


h=0 


d‘^T. 


s  d'^Tr 


dhidhj  dhidhj 


=  r  (4  -  4)  (4-81) 

^0  c^J  ■ 


Substituting  equations  (4.81),  (4.73)  and  (4.76)  into  (4.71),  and  using 


I  obtain 


dhidhj 


h=i 


2z 


d^z 


dhidhj  J 


h=0 


dhidhj 


h=k 


-2  I  2 


(4.82) 


(4.83) 


Finally,  by  using  the  above  equation,  the  residual  moveout  (4.69)  is  simplified  to 


z‘^{h)  =  z‘^{0)  +  -  ij  |/i|2  +  0(|/i|^),  (4.84) 

or 

^m(^)  =  rl{0)  +  4  Q  -  ^)  |/i|2  +  0(|/i|"),  (4.85) 

where  is  the  migration  time.  Equation  (4.85)  holds  for  constant  velocities  and  an  arbi¬ 
trary  reflector,  and  shows  that  the  residual  term  in  the  the  migrated  data  is  independent 
of  both  reflector  dip  and  the  source-receiver  azimuth. 

When  the  medium  velocity  is  varied,  one  should  replace  v  by  although  this 

statement  is  precise.  Therefore,  equation  (4.85)  offers  a  theoretical  basis  for  estimating 
the  rms  velocity  from  the  3-D  residual  moveout. 


4.4  Limitation  on  a  Two-layer  Model 

Now,  I  design  a  test  to  show  the  limitation  of  a  velocity  estimation  that  uses  the 
hyperbolic  residual  moveout  (4.27).  The  model,  shown  in  Figure  4.2  consists  of  two 
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constant  velocity  layers,  Vi  and  V2,  and  two  flat  interfaces.  The  ray  starts  vertically  from 
the  midpoint  y  on  the  surface  and  intersects  the  flrst  interface  with  the  angle  ai  from 
the  normal;  then  transmits  to  the  second  layer  with  the  angle  /?i  from  the  normal  and 
intersects  vertically  the  second  interface  at  the  reflection  point  {x,z). 


Fig.  4.2.  Two-layer  model. 


The  ray  is  perpendicular  to  the  surface,  so  dtfdy  =  0.  This  fact  means  Vrmo  =  v^mo- 
From  Hubral  and  Krey’s  formula  (1980),  I  have 


=  'iP' 

''rmo  ^nmo 


1 


^ti  -1-  At2 


Uj  "t"  1^2  ^^2 


cos^  ai 

COS^  Pi 


(4.86) 


where  Ati  and  A<2  are  the  one-way  time  along  the  ray  in  the  flrst  and  the  second  layer 
respectively.  Suppose  that  is  10  percent  greater  than  Vi,  i.e.,  V2  =  l.lui.  Also,  I 
specify,  Ati  =  At2  and  Oi  =  60°.  After  calculations,  I  find  Pi  =  72.3°  and  Vrmo  =  1.46ui. 
From  the  Dix  equation,  the  estimated  interval  velocity  is  obtained. 


V2  =  1.81ui  =  1.64u2. 

That  is,  if  the  error  in  initial  velocity  is  10  percent,  then  the  error  in  the  updated  velocity 
may  be  magnified  to  64  percent. 

This  example  shows  the  limitation  of  velocity  analysis  by  hyperbolic  residual  move- 
out:  when  the  velocity  distribution  contains  a  strong  lateral  variation,  velocity  updating 
may  fail  to  improve  the  initial  value. 

4.5  Computer  Implementation 

Parallel  to  NMO  velocity  analysis,  one  can  use  a  semblance  and  velocity  scans  to  do 
RMO  velocity  analysis.  Therefore,  the  program  for  RMO  velocity  analysis  is  similar  to 
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the  existing  program  for  NMO  velocity  analysis,  except  in  the  residual  term.  A  suggested 
data  processing  technique  based  on  this  approach  is  composed  of  prestack  migration  with 
a  constant  velocity,  velocity  analysis,  residual  moveout,  stacking,  and  poststack  residual 
migration. 

To  demonstrate  the  processing  procedure,  I  applied  this  method  to  synthetic  seismo¬ 
grams  computed  for  a  subsurface  model  in  which  velocity  increases  linearly  with  depth 
according  to  v{z)  =  1.5-1-  O.82  km/s.  The  model,  shown  in  Figure  4.3,  consists  of  five 
reflectors,  each  with  a  dipping  and  horizontal  segment.  Dips  for  the  dipping  segments 
range  from  30  to  90  degrees  in  15-degree  increments.  The  seismograms  contain  10  ofi"- 
sets,  ranging  from  100  to  1900  meters  in  200  meter  increments.  Because  of  the  dipping 
reflectors  and  the  depth  dependent  velocity,  the  RMO  velocity  in  equation  (4.53)  and 
the  rms  velocity  are  not  the  same  but  close  each  other.  The  relative  error  between  the 
two,  shown  in  Figure  4.4,  increases  with  depth  and  dip.  The  maximal  error  is  about  two 
percent. 

After  prestack  migration  with  the  constant  velocity,  c  =  1.5  km/s,  the  dipping  events 
are  not  migrated  to  correct  positions,  except  for  the  30-degree  dip  (see  Figure  4.5).  One 
of  the  common  image  gathers  is  plotted  in  Figure  4.6.  Because  the  migration  velocity 
is  lower  than  the  true  velocity,  the  imaged  depths  increase  with  offset.  The  velocity 
spectra  for  this  CIG  is  shown  in  Figure  4.7.  Unlike  the  velocity  spectrums  in  NMO, 
the  velocity  peaks  here  are  insensitive  to  reflector  dips  and,  therefore,  single- valued. 
The  RMO  velocities  are  picked  from  these  peaks.  After  residual  moveout  correction, 
all  events  in  the  CIG  are  corrected  to  horizontal  ones,  shown  in  Figure  4.8.  Stacking 
the  data,  having  residual  moveout  corrected,  yields  a  result  shown  in  Figure  4.9,  which 
is  equivalent  to  the  zero-offset  migration  with  the  constant  velocity  c.  By  using  the 
interval  velocity  converted  from  the  RMO  velocity,  poststack  residual  migration  gives 
the  corrected  reflector  positions  shown  in  Figure  4.10.  Notice  that  the  bottom  of  the 
vertical  event  is  not  migrated  completely,  which  is  due  to  some  errors  involved  in  the 
RMO  velocity  estimates  and  that  these  errors  are  enlarged  in  the  conversion  of  the 
interval  velocity.  For  example,  the  relative  error  between  the  true  interval  velocity  at  2 
seconds  and  the  computed  one  is  4  percent. 

4.6  Summary 

Under  the  assumption  of  small  offset,  residual  moveout  representations  are  obtained 
in  this  chapter.  Using  these  relationships,  one  can  estimate  directly  the  medium  velocities 
from  residual  moveouts  without  iteration  when  the  medium  has  weak  lateral  velocity 
variations.  For  complex  structures,  this  approach  does  not  offer  an  accurate  enough 
velocity  estimation,  even  with  iteration. 

For  converted  waves  or  anisotropic  media,  seismic  data  are  not  symmetric  to  offset 
generally.  Thus,  the  residual  moveout  representations  will  be  dip-dependent  and,  there¬ 
fore,  be  less  useful  in  velocity  estimation.  However,  these  representations  may  allow  us  to 
make  a  qualitative  analysis  of  the  relationship  between  the  imaged  depth  and  migration 
velocity. 
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Fig.  4.3.  Subsurface  model  used  to  generate  synthetic  seismic  traces. 


Dip  (degree) 


Fig.  4.4.  The  contours  of  relative  error  between  the  RMO  velocity  and  the  rms 
velocity.  The  error  increases  with  depth  and  dip. 
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4.5.  The  migrated  data  with  the  contant  velocity.  The  offset  is  100  meters. 


Offset  (km) 
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Fig.  4.6.  One  of  the  common  image  gathers  from  date  migrated  with  the  constant 
velocity.  The  traces  are  located  at  1.4  km  in  Figure  4.5. 
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Fig.  4.7.  Velocity  spectrum  for  the  CIG  in  Figure  4.6. 
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Fig.  4.8.  Residual  moveout  correction  for  the  CIG  in  Figure  4.6. 
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Fig.  4.9.  Stacking  for  the  partially  migrated  data  of  the  ten  offsets. 


Midpoint  (km) 


Fig.  4.10.  Poststack  residual  time  migration  of  the  data  in  Figure  4.7. 
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Chapter  5 


ITERATIVE  APPROACH:  METHODOLOGY 


When  velocity  has  lateral  variations,  the  residual  moveout  cannot  be  approximated 
by  a  hyperbola  and  the  RMO  velocity  may  be  quite  different  from  the  rms  velocity. 
Therefore,  the  velocity  estimate  cannot  be  simply  done  by  using  residual  moveout  cor¬ 
rection.  In  this  situation,  iterative  approaches  are  required  to  update  velocity.  However, 
iterative  formulas  used  in  conventional  approaches  are  derived  under  assumptions  such 
as  small  offset,  small  dip,  and  lateral  velocity  homogeneity.  Although  iteration  generally 
is  helpful  in  obtaining  a  more  accurate  velocity,  too  coarse  an  approximate  formula  for 
updating  velocity  not  only  increases  the  number  of  iteration  steps  but  may  result  in  di¬ 
vergence.  Thus,  applications  of  those  approaches  to  velocity  analysis  are  limited  when 
there  are  complex  structures.  In  Chapter  4, 1  have  shown  a  counter  example  for  iteration 
when  the  hyperbolic  residual  moveout  method  is  applied  to  a  simple  two-layer  model. 

In  this  chapter,  I  will  present  an  iterative  approach  to  migration  velocity  analysis. 
Using  perturbation  theory,  I  derive  a  formula  for  updating  velocity  from  residual  moveout 
by  computing  a  derivative  function  of  imaged  depths  with  respect  to  velocity.  Signifi¬ 
cantly  here,  this  formula  has  no  limitations  on  offset,  reflector  dip,  or  velocity  distribution 
if  the  velocity  perturbation  is  sufficiently  small.  Therefore,  this  formula  gives  a  general 
description  of  the  relationship  between  residual  moveout  and  residual  velocity.  Based  on 
this  formula,  I  revise  the  residual-curvature- analysis  method  for  velocity  estimation.  In 
addition,  this  formula  provides  both  sensitivity  and  error  estimation  for  migration-based 
velocity  analysis,  which  is  helpful  in  explaining  the  reliability  of  the  estimated  veloc¬ 
ity.  Moreover,  this  velocity  analysis  approach  can  be  extended  to  the  3-D  case  and  to 
converted  waves. 

5.1  Velocity  Analysis  by  Perturbation 

In  Chapter  4,  the  residual  moveout  representations  are  obtained  under  the  small- 
offset  assumption.  However,  this  assumption  is  violated  for  a  lateral  velocity  variation. 
Here,  I  will  derive  a  residual  moveout  representation  by  using  a  perturbation.  This  idea 
is  suggested  by  tomography  approaches  to  velocity  estimation.  Under  the  assumption  of 
small  velocity  perturbation,  there  is  a  linear  relationship  between  the  residual  traveltime 
and  the  residual  velocity.  In  the  tomographic  approach,  raypaths  are  required  to  con¬ 
struct  the  equation  system  for  velocity  estimation.  However,  raypath  construction  seems 
more  difficult  in  reflection  tomography  than  in  refraction  tomography.  Raypaths  strongly 
depend  on  the  reflector  position  (especially  the  slope)  in  reflection  tomography  and  the  re¬ 
flector  cannot  be  determined  accurately.  In  contrast,  here  I  will  use  the  stationary-phase 
principle  to  determine  raypaths,  without  requiring  knowledge  of  an  accurate  reflector 
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position.  My  derivation  will  be  based  on  the  imaging  equations  in  Chapter  2. 


5.1.1  Mathematical  derivation 

Consider  the  2-D  common  offset  case.  All  notations  have  the  same  meaning  as  those 
in  Chapter  2.  The  imaging  equations  (2.8)  and  (2.9)  are  repeated 

T^{X^,X)  +  Tr(x,Xr)  = 

OTs  dXr  _  dt 

dy  dy  dy' 

At  a  common  image  gather,  the  imaged  depth  z  can  be  determined  as  a  function  of  h.  If 
the  migration  velocity  equals  the  true  velocity,  then  the  imaged  depth  z  will  be  indepen¬ 
dent  of  offset  h\  otherwise — for  incorrect  velocity — z  varies  with  offset  h.  Consequently, 
the  imaged  depths  in  CIGs  provide  information  on  velocity  distribution. 

Equations  (5.1)  and  (5.2)  display  a  general  relationship  between  the  imaged  depth 
and  migration  velocity.  However,  this  equation  system  is  nonlinear,  making  it  very 
difficult  to  directly  solve  for  velocity.  Here,  I  use  a  mathematical  tool — perturbation — to 
linearize  this  equation  system. 

Suppose  that  the  velocity  distribution  v  is  characterized  by  a  parameter  or  a  family 
of  parameters.  A, 

V  =  v{x\  A). 

For  example,  when  v{x',  A)  =  uq  +  oa;  -f  fez,  A  is  any  set  of  one  to  three  parameters 
chosen  from  vq,  a,  and  fe.  Thus,  the  problem  of  velocity  estimation  becomes  parameter 
estimation.  To  simplify  the  derivation,  I  suppose  that  A  is  just  a  single  parameter  at 
first. 

For  a  fixed  image  location  x,  I  differentiate  equation  (5.1)  with  respect  to  A.  Noticing 
that  y  and  z  are  functions  of  A;  then 


(5.1) 


dxg  dxr  dy  dxg  dr  A  \dTa  dz  _  dt  dy 

dy^  dy  dX"^  d\^  dX/^  dz^  dz  dX  dydX’ 


By  using  equation  (5.2),  the  first  term  of  the  left  side  in  equation  (5.3)  is  balanced  by 
the  right-hand  term.  Therefore, 


Let  6s  or  6^  be  the  angle  between  the  raypath  from  the  source  or  the  receiver,  and  the 
vertical  at  x\  then 

dxs  cos 6s  dxr  cos6r 


v(x;X)' 


v(x;X)’ 
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By  using  the  above  equation,  equation  (5.4)  is  rewritten  as 


cos  9a  +  cos  6a  dz 


dVa  dTr 

'  dX  dX' 


v{x',X)  dX 

Thus  the  derivative  of  imaged  depth  with  respect  to  A  is  found, 

dz 


dX 


=  g{x,h), 


(5.6) 


(5.7) 


where 


g{x,h)  =  - 


dta  dtr 

'dX'^Jx 


v{x]X) 

cos  9a  +  cos  9r 


(5.8) 


The  function  g  characterizes  the  relationship  between  the  imaged  depth  and  the  migration 
velocity  in  a  general  medium  context.  The  computation  of  this  function  will  result 
in  a  new  migration  velocity  analysis  method,  compared  to  conventional  ones  based  on 
hyperbolic  residual  moveout. 

Suppose  that  the  true  parameter  is  A*  and  the  true  reflection  depth  is  z*.  If  there 
is  a  small  perturbation  6X  =  X*  —  X  between  the  true  parameter  and  the  parameter  used 
in  migration,  then  the  imaged  depth  will  have  a  corresponding  perturbation 


6z{x,  h)  =  z*  —  z{x,  h) 


By  using  equation  (5.7),  a  linearized  equation  is  set  for  8X\ 


dz  =  g{x^  h)6X. 


(5.9) 


Equation  (5.9)  is  valid  for  an  arbitrary  velocity  distribution,  arbitrary  reflector  dips, 
and  any  offset,  as  long  as  the  velocity  perturbation  is  sufficiently  small,  which  is  signifi¬ 
cantly  different  from  the  limited  result  of  conventional  RCA. 

Remark.  When  the  velocity  distribution  is  constant  in  layers,  one  can  take  A  as 
velocity  v  itself,  in  the  target  layer.  In  this  situation. 


dXa  dTr  _  ta+tr 

dX^  dX~  v 


(5.10) 


where  tg  (or  tr)  is  a  partial  traveltime  of  (or  r^)  within  this  layer.  Thus,  the  derivative 
function  g  is  simplifled  to 


g{x,h) 


tg  ”1"  tr 

COS  63  +  cos  9r 


(5.11) 


This  typical  result  was  obtained  by  Lafond  and  Levander  (1993). 

When  the  velocity  distribution  is  characterized  by  multiple  parameters,  A  =  (Ai,  A2, 
A„)^,  the  imaged-depth  perturbation  will  depend  on  the  perturbations  of  all  these 
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parameters.  Therefore,  equation  (5.9)  will  be  modified  by 

"  8z  ” 

6z{x,h)  =  =  Y^gi{x,h)SXi, 

i=l  ,=1 

where  each  derivative  function  is  calculated  by 


dXi  dXi  cos  +  cos  ffr 


(5.12) 


(5.13) 


If  one  could  solve  for  5A,  the  true  parameters  can  be  estimated  by 

y  =  X  +  SX.  (5.14) 

Notice  that  the  left  side  of  equation  (5.12)  involves  the  true  depth  z*  that  is  unknown. 
Conventional  methods  deal  with  this  problem  by  using  the  concept  of  the  reference  depth 
(Lafond  and  Levander,  1993)  or  adding  z*  as  a  new  unknown  (Stork,  1991).  Here,  I  use  a 
technique  to  remove  z*  directly  in  equation  (5.12).  The  true  depth  can  be  approximately 
replaced  by  the  corrected  imaged  depth  z  +  Sz, 

n 

z*  «  z(x,  h)  +  6z{x,  h)  =  z{x,  h)  +  gi{x,  h)6Xi.  (5.15) 

i=l 

The  true  reflection  depth  z*,  is  independent  of  offset.  Therefore,  the  corrected  imaged- 
depths  from  different  offsets  should  be  close  to  each  other.  Mathematically,  this  statement 
means  that  the  covariance  of  these  depths  is  a  minimum.  Suppose  that  there  are  offsets 
hi,  /i2)  •••)  hm,  and  image  locations  Xi,  X2,  ...,  xk,  then 

zf  +  6z^'‘^  =  zf  +  E4^^<5A.-,  (5.16) 

i=l 

where 

zf^  =  z{xk,hj), 

Sz^^^  =  5z{xk,hj), 

9ij^  =  9i{xk,hj). 

I  seek  (5A,’s  such  that  the  corrected  imaged  depths  have  the  minimum  variance;  i.e., 

K  m  _  ty 

{zf^  +  (5zj^^  —  z(^)  +  6z(*) j  =  min,  (5-17) 

ifc=ii=i  ^ 


where 


z(*^)  =  (zf\zf,...,zW)^ 
6z^^^  =  {6zt\6zP,...,6z^,^^f. 
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Here  I  use  the  overline  to  denote  the  mean  value  of  a  vector  over  the  offset  index.  For 
example, 


Indroduce  the  matrix  and  vector, 


where 


(5.18) 

where 

i=i  ^  ^ 

.(k)  /  (fc)  (k)  (k)xT 

9i  =  {gll  ,9i2  ,-^9in  ) 

then,  as  shown  in  Appendix  A,  the  solution  of  equation  (5.17)  must  satisfy  the  linear 
equation, 


'  AT  1  iC 

Jc=l  J  fc=l 


(5.19) 


Specifically,  if  there  is  only  one  parameter  to  be  determined,  i.e  n  =  1,  then  equa¬ 
tion  (5.19)  will  have  an  explicit  solution 


where 


Ef.i  Z",!  (gf  -  (#’  - 

ELi  E” ,  (sf  - 1®)" 


9j  =  9(^k,hj), 


(5.20) 


If  the  corrected  imaged-depths  are  not  close  enough  to  each  other,  I  implement  iter¬ 
ation  to  obtain  more  accurate  parameters.  The  iteration  stops  where  the  variance  (5.17) 
achieves  a  given  sufficiently  small  value. 


5.1.2  Computational  techniques 

Calulation  of  the  derivative  function 

The  function  g(x,  h)  involves  the  derivatives  of  traveltimes  with  respect  to  the  pa¬ 
rameter  A.  In  the  eikonal  equation. 


dx)  ^  \dzj 


(5.21) 


DiflFerentiating  here  with  respect  A  yields 


dfidr  dfidr  _  1  d  /  1  \ 

dx  dx  dz  dz  2  dX  y  A)  J  ’ 


where 


The  integral  solution  of  equation  (5.22)  is 


H  = 


(5.22) 


(5.23) 


where  L  is  the  raypath  from  the  source  (or  receiver)  to  the  image  point  x. 

For  each  source  or  receiver,  drfdX  can  be  determined  from  equation  (5.22)  or  (5.23). 
Therefore,  given  an  image  point  x  and  a  specular  source-receiver  pair  Xg  and  Xr,  one  can 
calculate  g  from  formula  (5.8).  However,  there  is  not  an  explicit  formula  to  represent  the 
specular  source-receiver  pair  from  the  image  point  for  a  complex  medium.  To  solve  this 
problem,  here,  I  use  the  Kirchhoff  integral  to  calculate  g.  In  the  Kirchhoff  summation, 
I  calculate  two  migration  outputs  which  have  the  same  phase  but  different  amplitudes. 
One  uses  the  original  amplitude;  the  other  one  uses  the  original  amplitude  multiplied 
by  the  quantity  g.  Thus,  the  ratio  of  the  amplitudes  of  these  two  outputs  will  evaluate 
g  at  the  specular  source-receiver  position  according  to  the  stationary-phase  principle, 
without  requiring  knowledge  of  the  specular  source-receiver  pair.  This  technique  is  the 
same  as  was  used  to  determine  the  angle  of  reflection  in  Kirchhoff  inversion  (Bleistein  et 
al.,  1987). 

This  approach  to  calculate  the  derivative  function  g  depends  the  S/N  ratio  in  seismic 
data  and  the  power  of  stationary  phase  method.  The  stationary  phase  method  works 
well  if  the  dominant  seismic  wavenumber  is  large  compared  to  the  length  scale  of  the 
velocity  variation.  In  this  sense,  a  smooth  velocity  is  required  in  Kirchhoff  migration  for 
velocity  estimates  by  perturbation. 


Parameterization  of  velocity  distribution 

Although  equation  (5.19)  holds  for  any  velocity  distribution,  the  solution  will  be  un¬ 
derdetermined  and  unstable  if  too  many  unknown  parameters  are  involved.  Consequently, 
it  is  essential  to  characterize  the  velocity  distribution  by  choosing  appropriate  parame¬ 
ters.  Conventionally,  one  assumes  that  a  velocity  model  consists  of  the  construction  of 
the  macro-model  (constant  velocities  and  velocity  interfaces) .  The  interfaces  divide  the 
whole  model  into  a  number  of  blocks  (shown  in  Figure  5.1).  Here,  I  replace  constant 
velocity  in  one  block  by  a  linear  function  that  is  characterized  by  three  parameters: 


Al  +  \2{z  -  ^o)  +  Xz{x  -  Xo), 
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where  {xo,Zo)  is  a  reference  point.  Thus,  the  velocity  distribution  is  written  in  a  form 
v{x,  z)  =  uo(a:,  z)  +  Ai  +  X2{z  —  zq)  +  X3{x  —  xq),  (5.24) 

where  uq  is  a  background  velocity. 

An  imaged  depth  only  depends  on  the  velocity  above  it,  except  for  turning  rays. 
Therefore,  use  of  a  recursive  algorithm  (layer  stripping)  is  possible  to  determine  ve¬ 
locity  in  an  individual  block.  I  start  from  the  block  nearest  surface.  For  example,  1 
choose  the  top  left  block  in  Figure  5.1.  In  each  block,  iteration  is  used  to  calculate 
velocity  parameters.  Given  an  initial  guess  for  A,’s,  common-offset  depth  migration  is 
implemented  to  obtain  imaged  depths  and  ffi(x,h)  in  equation  (5.13)  for  common  image 
gathers.  Equation  (5.19)  will  give  a  correction  of  the  parameters.  Then  by  using  the 
updated  parameters  as  an  initial  guess,  I  correct  the  velocity  again  until  convergence  is 
achieved.  After  velocity  analysis  in  one  block,  I  migrate  data  with  the  corrected  velocity, 
and  pick  the  velocity  interface  from  the  imaged  structure.  When  I  finish  determining 
velocity  and  velocity  interface  in  one  block,  I  will  repeat  the  same  procedure  to  the  next 
block  that  is  located  below  the  finished  block. 

An  alternative  choice  of  perturbation  approach 

Perturbation  methods  are  used  to  gain  good  linear  approximations  to  nonlinear 
functions.  However,  in  some  cases,  one  should  exercise  caution  in  choosing  the  variable 
in  which  an  expression  is  perturbed.  For  example,  in  acoustic  media,  z^  as  a  linear 
function  of  the  square  of  migration  velocity  is  more  likely  than  z  as  a  linear  function 
of  the  migration  velocity  itself.  For  example,  for  a  constant  velocity  medium  and  a 
horizontal  reflector,  equation  (4.3)  shows  that  ^  is  the  exact  linear  function  of  the  square 


45 


(5.25) 

(5.26) 

(5.27) 


5.2  Sensitivity  of  Migration  Velocity  Analysis 

Velocity  analysis  by  prestack  migration  uses  the  difference  between  the  imaged 
depths  from  different  offsets  to  correct  the  velocities,  which  is  represented  by  equa¬ 
tion  (5.19).  If  the  variance  defined  by  equation  (5.17)  is  zero,  one  can  conclude  that 
the  velocity  is  correct.  However,  it  is  impractical  to  obtain  an  exact  zero  variance.  There 
are  many  factors  against  achieving  this  goal:  noise  in  the  input  data,  nonacoustic  prop¬ 
erties,  inaccurate  description  of  velocity  distribution,  and  so  on.  Even  while  an  apparent 
zero  variance  is  obtained,  the  actual  one  may  not  be  zero  because  of  the  errors  involved 
in  picking  the  imaged  depths.  The  imaged  depths  in  the  variance  are  picked  on  the  mi¬ 
gration  output,  so  that  the  position  error  in  imaged  depths  is  controlled  by  the  resolution 
of  the  migration  output,  which,  in  turn,  depends  on  wavelength,  among  other  things. 

Quantitatively,  the  matrices  in  equation  (5.19)  that  depend  on  the  functions  p,,  can 
be  used  to  describe  the  sensitivity  of  the  velocity  error  ^A,  to  the  variance  of  residual- 
moveout  error.  Here,  I  will  derive  analytical  representations  for  the  simplest  cases. 

Suppose  that  the  velocity  v{x,  z)  consists  of  a  constant  background  velocity  vq  and 
a  perturbation  that  is  a  linear  function  of  depth  in  one  block.  Moreover,  I  assume  that 
the  upper  boundary  of  the  block  is  a  horizontal  line,  z  =  d,  and  the  reflector,  located  at 
z  =  z*,  is  the  bottom  boundary  of  this  block.  The  velocity  function  can  be  represented 
by 

v{x,  z)  =  Vo  +  q;(Ai  +  \2{z  -  Zq)),  (5.28) 

where  a  =  0  for  z  <  d,  a  =  1  for  z  >  d,  and  zq  is  a  reference  depth  (not  necessarily  in 
the  block).  A  sketch  for  this  block  is  shown  in  Figure  5.2. 

The  initial  guesses  are  set  to 

Ai  =  A2  =  0. 


As  shown  in  Appendix  B,  I  obtain  the  following  representations  for  the  p,  ’s  that  is  defined 


in  equation  (5.13), 


gi{x,h) 


+  z"^  z  —  d 


(5.29) 
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Fig.  5.2.  Layer  model.  The  gray  shading  denotes  the  target  layer  for  velocity  analysis. 


g2{x,h) 


h?  +  z‘^  (z  —  zo)^  —  (d  —  zq)^ 
2vo 


From  the  above  two  equations,  I  have 


92{x,h) 


{z  -  Zpf  -{d-  Zpf 
2{z  —  d) 


gi{x,h). 


(5.30) 


(5.31) 


This  means  that  the  coefficients  in  equation  (5.12)  are  proportional  for  all  image  locations 
and  all  offsets,  so  that  I  cannot  solve  for  Ai  and  A2  independently.  Thus,  I  have  the 
following  conclusion:  The  parameters  Ai  and  A2  cannot  he  determined  simultaneously  if 
only  one  reflector  segment  is  used  in  the  velocity-analysis  process. 

This  conclusion  shows  that  one  should  avoid  solving  for  Ai  and  A2  at  the  same  time 
unless  well-separated  reflector  segments  are  used  for  velocity  analysis. 

If  the  parameter  A2  is  given  (i.e.,  6X2  =  0),  then  equation  (5.13)  is  simplified  to 


(5.32) 


When  (5Ai  is  small,  z{x,h)  «  2*.  Therefore,  for  given  two  offsets  ^2  and  hi  (^2  >  hi), 
the  difference  of  imaged  depths  from  these  two  offsets  is 


z{x,h2)  -  z{x,hi) 


^1 

(Z*)2_ 


z*  —  d 


vq 


6\i. 


(5.33) 
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(5.34) 


I\Z  —  Z\X^h2)  (yAj. 

\z^rV{i 

The  above  equation  shows  that  the  error  in  Ai  is  proportional  to  the  thickness  of  the  block 
layer,  which  is  consistent  with  results  of  equation  (3.34).  When  there  is  a  thin  layer,  the 
velocity  estimation  for  this  layer  is  often  unreliable,  so  I  conclude  that  migration-based 
velocity  analysis  cannot  handle  thin  layers. 

Similarly,  if  the  parameter  Aj  is  given  (i.e.,  <5Ai  =  0),  then 


Az  = 


{h^  —  hi)(2*  —  d){z*  +  d  —  2zo) 
2{z*)^vo 


6X2. 


(5.35) 


The  above  equation  shows  that  the  error  in  A2  is  determined  not  only  by  the  thickness 
of  the  block  layer  but  by  {z*  +  d)f2  —  Zq — the  difference  between  the  central  depth  of  the 
block  layer  and  the  reference  depth  Zq. 

There  are  two  kinds  of  depth  errors:  6z  and  Az.  6z  is  the  difference  between  the  true 
depth  and  the  imaged  depth,  which  reflects  migration  error  due  to  the  velocity  error;  Az 
is  the  difference  between  the  imaged  depths  from  different  offsets,  i.e.,  residual  moveout. 
The  ratio  of  these  two  depth  errors  is  given  by  equations  (5.32)  and  (5.34), 

^  +  (z*f 

^  ~  Az  (/i2  —  hi) 

The  ratio  7  indicates  the  depth  error  of  migration  due  to  the  error  in  migration  residual 
moveout.  The  bigger  the  value  of  7  is,  the  larger  error  of  migration  may  result  from 
the  error  in  residual  moveout.  Equation  (5.36)  shows  that  increasing  the  ratio  of  offset 
to  reflection  depth  will  be  helpful  to  reduce  the  migration  error,  although  this  is  a  well- 
known  statement. 

For  a  general  case,  it  is  difficult  to  obtain  analytical  representations  for  the  ^,’s. 
However,  I  can  calculate  the  gfs  numerically,  as  described  above,  and  use  these  values  to 
estimate  the  velocity  error.  Specifically,  if  there  is  only  one  parameter,  equation  (5.20) 
gives  a  direct  error  estimation. 


(5.36) 


i«i< 


where  the  Cauchy-Schwarz  inequality  has  been  used. 


(5.37) 


5.3  Simple-Iteration  Approach 

The  velocity  inversion  essentially  is  a  nonlinear  process.  By  using  the  perturbation 
method,  I  obtain  a  linear  equation  (5.19)  to  update  velocity  iteratively.  The  efficiency  of 
this  linearized  iteration  depends  on  the  initial  guess  of  velocity  and  the  stability  of  velocity 
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updating.  Actually,  these  problems,  initial  guess  and  stability,  are  related.  The  less  the 
stability  is,  the  better  the  initial  guess  that  is  required.  The  stability  of  equation  (5.21) 
depends  on  the  noise  level  in  the  data  and  the  complexity  of  structures  in  the  overburden. 
For  example,  equation  (5.21)  trends  unstable  when  a  thin  layer  structure  exists.  In  order 
to  obtain  a  convergent  solution,  I  use  the  so-called  simple-iteration  algorithm  decribed 
below  as  a  supplement  to  the  perturbation  method.  Here,  I  assume  that  the  velocity  is 
constant  layered. 

If  I  view  residual  moveout  Az,  as  a  function  of  migration  velocity  v,  then  the  true 
velocity  v*  is  the  one  for  which 

Az{v)  =  0.  (5.38) 

Also,  from  equation  (5.34),  I  conclude  that  Az  is  positive  when  v  >  v*  and  Az  is 
negative  when  v  <  v*.  Based  on  this  fact,  one  proceeds  with  the  following  simple- 
iteration  algorithm  to  determine  the  true  layer  velocity  as  follows: 

(1)  Choose  two  initial  velocities  v\  and  V2  such  that 

Az{v\)  <  0,  Az{v2)  >  0. 

(2)  Compute  a  new  velocity  by  weighting  the  initial  velocities  as  follows: 

—  Az{v2)  Az{vi) 

^  Az{V2)  —  Az{vi)  Az{vx)  —  Az(v2)  ' 

(3)  If  Az{v3)  =  0  (or  smaller  than  a  given  precision),  the  iteration  will  stop,  with  V3 
the  desired  velocity.  If  Az{v3)  >  0,  V3  replaces  V2  in  step  (1);  otherwise,  if  Az{v3)  <  0, 
V3  replaces  vi  and  go  to  step  (2). 

During  the  iteration,  the  true  velocity  v*  always  is  between  vi  and  V2,  and  the 
deviation  of  v\  and  V2  decreases  monotonically,  so  the  convergence  of  this  iteration  is 
guaranteed  as  long  as  the  model  is  well  approximated  by  constant-velocity  layers. 

5.4  Some  Extensions 

It  is  straightforward  that  formulas  in  section  4.1  are  extended  from  the  common 
offset  case  to  the  common  shot  case,  because  I  did  not  use  the  special  properties  of 
common  offset  gather  in  my  derivation.  The  only  difference  is,  traces  in  a  CIG  for  the 
common  shot  case  are  organized  in  accordance  with  shot  positions  instead  of  the  source- 
receiver  offsets.  Now  I  will  discuss  the  extension  of  these  results  to  the  3-D  case  and 
converted  waves. 

5.4.1  3-D  case 

The  derivation  in  3-D  is  the  same  as  that  in  section  4.1.1,  except  for  replacing 
location  scalars  by  location  vectors  and  replacing  offset  scalars  by  offset  vectors.  For 
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example,  equation  (5.8)  will  be  rewritten  as 


(5.40) 


g{x,h)  =  - 


dr,  dtr 

dX  dX 


vjx^X) 

cos  9s  +  cos  Or  ’ 


Here  x  =  {xi,X2),  and  h  =  {hi,h2)  have  the  same  meaning  as  that  in  section  2.3.  In 
addition,  the  linear  velocity  function  will  depend  on  four  parameters: 


v{x,  z)  =  uo(ar,  2)  +  Ai  +  A2(^r  -  zq)  +  A3(a:i  -  xio)  +  A4(a;2  -  a;2o),  (5.41) 

where  vq  is  a  background  velocity  and  (xio,  X20,  ^o)  is  a  reference  point. 

Theoretically,  the  velocity  analysis  approach  in  this  chapter  will  be  suitable  both 
to  the  2-D  and  to  the  3-D  case.  However,  this  approach  certainly  encounters  more 
computational  difficulties  in  the  3-D  case.  Two  issues  are  (i)  how  one  measures  residual 
moveout  efficiently,  and  (ii)  how  one  picks  the  imaged  interface  in  the  migrated  data  for 
building  a  macro  model.  Hopefully,  these  difficulties  will  be  solved  with  development  of 
3-D  data  processing  and  visualization. 


5.4.2  Converted  waves 

With  some  work,  the  formulas  in  Section  4.1  can  be  extended  to  converted  waves.  In 
each  expression,  there  is  no  difference  between  imaging  equations  for  converted  waves  and 
those  for  non-converted  waves.  When  the  perturbation  method  is  applied  to  converted 
waves,  one  should  use  different  parameters  to  characterize  P-wave  velocity  vp,  and  S-wave 
velocity  U5,  even  within  a  same  block. 

With  no  loss  generality,  I  consider  a  P-S  refection  and  show  how  to  modify  some 
formulas  in  section  4.1.1.  Suppose  that  A  is  a  parameter  to  describe  then,  compared 
to  equation  (5.5),  I  have 


Bts  cos  0. 


dtr  cos  6r 


dz  vp{x-,Xy  dz  vs{x‘,X)‘ 

Therefore,  compared  to  (5.8),  the  function  g  has  a  modified  representation. 


g{x,h)  = 


dxr 


cos  9.  cos  9r 


n-i 


+ 


dX  [up(a;;A)  us(aj;A)J 


Here  I  use  the  fact  that 


Parallel  to  equation  (5.23),  I  have 

dxr 

lx 


S-” 


■/-I— 

Jl  dX  \vs{x; 


dL, 


(5.42) 


(5.43) 


where  L  is  the  raypath  from  the  receiver  to  the  image  point  x. 
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(5.44) 


5.5  Summary 


Perturbation  methods  are  general  approaches  to  solving  nonlinear  mathematical 
problems.  Based  on  perturbation  methods  and  the  imaging  equations,  I  derived  a  rela¬ 
tionship  between  residual  moveout  and  residual  velocity  that  has  no  limitations  on  reflec¬ 
tor  dip,  offset  or  velocity  distribution.  The  proposed  algorithm  estimates  the  update  in 
velocity  by  computing  a  derivative  function  of  imaged  depths  with  respect  to  velocity  in 
a  general  background  medium  context.  This  formula  is  more  accurate  than  conventional 
formulas  based  on  hyperbolic  residual  moveout  when  the  medium  has  strongly  lateral  ve¬ 
locity  variations.  From  this  point  I  have  revised  the  conventional  RCA  approaches.  With 
a  proper  modiflcation,  the  formulation  here  is  also  suitable  to  the  3-D  case,  converted 
waves,  and  even  anisotropic  media.  However,  for  anisotropic  media,  one  requires  more 
parameters  to  describe  a  velocity  distribution,  so  the  equation  corresponding  to  (5.19) 
will  be  underdetermined  if  only  one  component  of  residual  moveout  is  used.  In  any 
case,  the  methodology  here  provides  a  basis  for  developing  computational  techniques  of 
velocity  analysis  in  anisotropic  media. 
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Chapter  6 


ITERATIVE  APPROACH:  COMPUTER  IMPLEMENTATION 


To  demonstrate  the  effectiveness  and  the  efficiency  of  the  velocity  analysis  techniques 
in  Chapter  5,  I  present  some  numerical  examples  for  common-offset  experiments.  The 
experiments  include  synthetic  data,  physical-tank  data,  and  Marmousi  data.  For  these 
experiments,  I  use  the  Kirchhoff  integral  to  implement  prestack  depth  migration  (Liu, 
1993).  The  layer-stripping  procedure  for  velocity  analysis  can  be  stated  as  follows: 

•  begin  from  the  first  block 

1.  estimate  velocity  parameters  iteratively 

(a)  migrate  with  an  initial  guess  of  velocity; 

(b)  sort  the  migrated  data  into  common  image  gathers; 

(c)  measure  imaged  depths  and  evaluate  the  derivative  function; 

(d)  update  velocity  by  using  the  perturbation  formula; 

2.  image  velocity  interface  by  using  corrected  velocity 

•  repeat  step  1  and  2  for  next  block 

When  velocity  and  velocity  gradients  are  estimated  simultaneously  in  one  block,  the 
iteration  tends  to  be  unstable,  as  explained  in  the  previous  chapter.  To  overcome  this 
difficulty,  it  is  preferable  to  estimate  velocity  first.  This  estimation  will  yield  an  averaged 
velocity  and  give  a  better  initial  guess  for  the  velocity  distribution  in  this  block. 

6.1  Synthetic  Data 

The  first  example  is  synthetic  data  generated  by  the  Kirchhoff  integral.  The  velocity 
model  shown  in  Figure  6.1  consists  of  six  blocks.  The  velocity  function  in  each  block  is 
constant  or  linear.  Synthetic  seismic  traces  are  generated  from  this  model,  with  five 
offsets  ranging  from  100  meters  to  900  meters.  Two  of  the  common  offset  gathers  are 
shown  in  Figure  6.2. 

The  velocity  analysis  process  is  outlined  as  follows: 

In  the  first  layer,  the  initial  guess  of  velocity  is  1500  m/s,  with  an  error  of  30  percent. 
Structural  images  are  distorted  in  the  migrated  data  using  this  erroneous  velocity,  shown 
in  Figure  6.3.  Two  common  image  gathers  of  the  migrated  data  also  indicate  that  the 
initial  guess  is  incorrect,  shown  in  Figure  6.4.  The  imaged  depths  are  measured  at  these 
two  image  locations  for  velocity  updating.  After  two  iterations,  the  corrected  velocity  is 
2004  m/s,  while  the  true  value  is  2000  m/s. 
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In  the  second  layer,  there  are  two  parameters  used  in  the  velocity  distribution:  a 
constant  velocity  and  a  vertical  constant  gradient.  The  initial  guess  takes  2004  m/s  that 
is  the  estimated  velocity  in  the  first  layer.  Four  common  image  gathers  of  the  migrated 
data,  shown  in  Figure  6.5  indicate  that  the  first  layer  velocity  in  the  initial  guess  is  correct, 
but  the  deeper  layer  velocities  are  not.  The  imaged  depths  are  measured  at  these  four 
image  locations  for  velocity  updating.  If  this  initial  velocity  is  directly  used  to  estimate 
the  velocity  constant  and  vertical  gradient,  the  iteration  will  be  divergent.  Therefore, 
I  only  estimate  the  velocity  constant  in  the  first  step.  The  updated  velocity  is  2420 
m/s  that  gives  a  better  approximation  to  the  true  velocity.  In  the  second  iteration,  the 
two  parameters  of  velocity  function  are  updated  simultaneously.  The  corrected  velocity 
function  is 

v(z)  =  2490  +  1.01(z  -  1000)  =  1480  +  l.Olz  m/s. 

The  true  velocity  function  is  v{z)  =  1500  +  l.Oz  m/s. 

In  the  left  third  layer,  there  are  three  parameters  used  in  velocity  distribution.  The 
initial  guess  is 

1480  +  1.012  =  2490  +  1.01(2  -  1000)  m/s, 

which  is  the  estimated  velocity  in  the  second  layer.  Three  common  image  locations  are 
used  for  updaing  velocity:  900,  1200  and  1500  m.  After  one  iteration,  the  corrected 
velocity  is 


v{x,  z)  =  2490  +  1.57(2  -  1000)  -  0.15x  =  920  +  1.572  -  0.15x  m/s. 

The  true  velocity  is  t;(x,  z)  =  1000  + 1.52  —  0.2a:  m/s.  Despite  the  error  in  the  lateral  gra¬ 
dient,  the  stopping  criterion  is  met  since  the  difference  of  the  imaged  depths  is  apparently 
zero.  Therefore,  iteration  ceases. 

In  the  right  third  layer,  there  are  two  parameters  used  in  velocity  distribution: 
constant  velocity  and  constant  vertical  gradient.  The  initial  guess  takes 

1480  +  1.012  =  2490  4-  1.01(2  -  1000)  m/s, 

which  is  the  estimated  velocity  in  the  second  layer.  Three  common  image  locations  are 
used  for  updating  velocity:  2550,  2850  and  3150  m.  After  one  iteration,  the  corrected 
velocity  is 

v{x,  z)  =  2490  1.55(2  -  1000)  =  940  +  1.552  m/s. 

The  true  velocity  is  ^(2)  =  1000  +  1.52  m/s. 

In  the  fourth  layer,  the  velocity  is  a  constant.  The  initial  guess  is  u  =  4000  m/s. 
Five  common  image  locations  are  used  for  updating  velocity:  1500,  1800,  1950,  2100  and 
2400  m.  After  two  iterations,  the  corrected  velocity  is  3534  m/s,  while  the  true  velocity 
is  V  =  3500  m/s. 

The  final  velocity  model  estimated  through  velocity  analysis  is  shown  in  Figure  6.6. 
Flat  residual  moveouts  at  common  image  gathers,  shown  in  Figure  6.7,  indicate  correct¬ 
ness  of  the  estimated  velocity.  With  this  velocity  model,  migration  is  implemented  and 
shown  in  Figure  6.8  which  is  close  to  the  migration  result  with  the  true  velocity,  shown 
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in  Figure  6.9. 


Table  6.1.  Test  for  ffi.  The  unit  of  offset  and  depth  is  m / s . 


Offset 

Imaged  depth 

Numerical  gx 

True  gx 

100 

370 

0.249 

0.251 

300 

360 

0.278 

0.282 

500 

330 

0.345 

0.346 

700 

290 

0.476 

0.475 

900 

210 

0.773 

0.783 

Also,  I  test  the  computation  of  gx  in  the  first  layer.  The  numerical  gx  is  computed 
from  the  ratio  of  the  amplitudes  in  the  Kirchholf  migration  outputs.  Then  I  compare 
the  numerical  gx  with  the  true  value  calculated  in  equation  (5.29).  In  this  test,  vq  =  500 
m/s  and  d  =  0.  The  result  listed  in  Table  (6.1)  shows  that  the  numerical  value  of  gx  is 
accurate  for  all  offsets.  This  validates  the  estimation  technique  based  on  the  stationary 
phase  principle. 

6.2  Physical-tank  Data 

The  input  data  is  from  a  physical  experiment  on  a  tank  model,  provided  to  us  by 
Marathon  Oil  Company.  The  real  medium  is  three-dimensional,  but  can  be  approximated 
to  a  two-and-half  dimensional  model.  Figure  6.10  shows  one  measured  slice  along  the 
in-line  direction.  The  data  were  296  shots,  each  shot  with  48  receivers.  The  shot  point 
spacing  is  24.38  m,  and  the  receiver  interval  is  24.38  m.  Five  common-offset  gathers  are 
used  for  velocity  analysis  whose  offsets  are  268.2,  512.1,  755.9,  999.7,  and  1243.6  meters. 
Two  of  them  are  shown  in  Figure  6.11.  The  first  shot  point  is  at  a:  =  0.  For  each  offset, 
there  are  256  shots  and  receivers.  The  time  sampling  interval  is  4  ms;  the  total  time  is  2  s. 
The  inversion  output  spans  the  ranges  x  from  61  to  6949  m  and  2:  from  0  to  3658  m.  In  this 
example,  the  “true”  velocities  do  not  exist,  since  the  model  is  not  perfectly  two-and-half 
dimensional.  Migration  velocities  are  sought  such  that  residual  moveout  is  minimized; 
i.e.,  migration  velocities  yield  the  best  stacking  effect  for  the  migrated  data.  The  property 
that  the  model  consists  of  constant-layered  velocities  allows  use  of  the  simple-iteration 
method  in  section  5.3  to  update  velocity.  Velocity  analysis  is  done  through  the  fifth  layer, 
and  the  first  four  ones  are  shown  in  Figure  6.12.  The  results  are  as  follows: 

(1)  In  the  first  layer,  the  iterative  values  of  cx  are  2438,  3962,  3309,  3570;  the 
measured  value  is  3581. 

(2)  In  the  second  layer,  the  iterative  values  of  C2  are  3570,  5182,  4779;  the  measured 
value  is  4801. 

(3)  In  the  third  layer,  the  iterative  values  of  C3  are  4779,  7315,  6681;  the  measured 
value  is  6831. 
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(4)  In  the  fourth  layer,  the  iterative  values  of  C4  are  6681,  3962,  4869;  the  measured 
value  is  4801. 

(5)  In  the  fifth  layer,  the  iterative  values  of  C5  are  4869,  9144,  5812,  6234;  the 
measured  value  is  6066. 

After  velocity  analysis,  a  velocity  model  is  obtained  that  consists  of  constant  veloc¬ 
ity  layers  and  interfaces,  shown  in  Figure  6.13.  The  stacked  migrated  data  using  this 
estimated  model  is  shown  in  Figure  6.14.  One  can  see  clearly  the  steep  dip  on  the  third 
interface,  the  saw-tooth  reflector  on  the  fourth  interface  and  the  flat  bottom  reflector  in 
this  figure.  In  contrast,  the  poorly  stacked  migration  section,  shown  in  Figure  6.15,  is 
obtained  with  the  measured  velocity  model  shown  in  Figure  6.10.  The  structural  images 
are  obviously  distorted.  Although  these  two  velocity  models  look  similar,  the  migration 
sections  show  significant  differences.  This  fact  implies  that  prestack  migration  is  very 
sensitive  to  the  velocity  model  and  migration  velocity  analysis  is  essential  for  imaging 
complex  structures. 

The  velocity  estimates  by  the  simple-iteration  use  only  one  common  image  gather. 
In  contrast,  the  perturbation  method  can  use  as  many  CIGs  as  possible  for  velocity  esti¬ 
mation.  Therefore,  the  perturbation  method  may  give  more  accurate  velocity  estimates. 
Here,  using  the  velocities  estimated  by  the  simple-iteration  as  an  initial  guess,  I  correct 
the  velocity  model  in  Figure  6.13  by  the  perturbation  method.  The  updated  velocity 
is  shown  Figure  6.16.  The  stacked  migrated  data  using  this  velocity  model,  shown  in 
Figure  6.17  gives  slightly  better  images  than  the  result  in  Figure  6.14. 

6.3  Marmousi  Data 

The  Marmousi  data  set  is  generated  by  using  a  two-dimensional  acoustic  finite- 
difference  modeling  program.  The  model  contains  many  reflectors,  steep  dips,  and  strong 
velocity  variations  both  in  lateral  and  vertical  directions  (with  a  minimum  velocity  of 
1500  m/s  and  a  maximum  velocity  of  5500  m/s),  shown  in  Figure  6.18.  The  data  set 
consists  of  240  shots  with  96  traces  per  shot.  The  initial  offset  is  200  m;  both  the  shot 
and  the  receiver  spacings  are  25  m.  The  first  shot  is  at  the  lateral  position  3000  m. 
Here  19  of  common-offset  data  gathers  are  used  for  velocity  analysis.  The  selected  offsets 
range  from  200  m  to  2000  m  with  spacing  100  m.  The  minimum-offset  gather  is  shown 
in  Figure  6.19. 

During  the  velocity  analysis  process,  I  assume  that  the  velocity  field  is  a  macro 
model  and  that  the  velocity  distribution  is  a  linear  function  in  each  block.  Velocity 
analysis  results  surely  depend  on  what  kind  of  migration  algorithm  is  used.  There  are 
two  commonly  used  approaches  to  calculate  traveltimes  in  Kirchhoff  migration:  finite 
differencing  and  ray  tracing.  Compared  to  ray  tracing,  the  finite  differencing  approach  is 
easier  to  code  and  more  efficient  to  implement,  but  fails  to  correctly  image  complicated 
structures  when  multiple  arrivals  exist  (Geoltrain  and  Brae,  1993). 

Here,  the  finite  differencing  is  initially  used  in  the  migration  implementation  for 
velocity  analysis  in  the  area  where  the  first  arrivals  carry  the  major  energy.  The  estimated 
velocity  model  is  shown  in  Figure  6.20.  In  the  central  bottom  parts,  a  paraxial  ray  tracing 
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algorithm  is  used  to  implement  Kirchhoff  migration.  In  this  approach,  the  traveltime 
corresponding  to  the  major  energy  is  chosen  when  multiple  arrivals  exist. 

Using  the  Kirchhoff  migration  algorithm  by  paraxial  ray  tracing,  velocity  analysis 
is  done  through  the  central  bottom  parts  of  the  Marmousi  model.  The  updated  velocity 
model,  shown  in  Figure  6.21,  consists  of  19  blocks.  In  each  block,  the  velocity  distribution 
is  a  constant  or  linear  function  of  the  depth.  Comparison  of  the  estimated  velocity  model 
to  the  true  one  at  three  lateral  locations  is  shown  Figure  6.22.  One  can  see  that  the 
estimated  velocity  matches  the  true  model  well  except  for  thin  layers.  In  fact,  from  the 
sensitivity  analysis  in  Chapter  5,  velocities  in  these  thin  layers  cannot  be  determined 
well.  The  stacked  migration  section  using  the  velocity  model  in  Figure  6.21  is  shown  in 
Figure  6.23.  Compared  to  the  migration  result  using  the  true  velocity  model,  shown  in 
Figure  6.24,  Figure  6.23  gives  an  acceptable  structural  image  even  in  the  central  bottom 
parts,  which  indicates  the  capability  of  this  migration  velocity  analysis  approach  for 
handling  complex  structures.  The  subsurfaces  are  well  imaged  except  for  some  detailed 
features  in  the  central  bottom  parts.  Some  blurry  image  in  the  central  bottom  parts 
may  be  caused  from  missing  high  velocity  zones.  Figure  6.18  shows  that  the  real  velocity 
contains  several  small  high-velocity  zones  in  the  central  parts.  These  high-velocity  zones 
do  not  appear  in  Figure  6.21  because  of  the  limitations  of  velocity  analysis  and  the 
resolution  of  migration  imaging. 

Selected  common  image  gathers  from  migrated  data  using  the  estimated  model  are 
shown  in  Figures  6.25,  6.26  and  6.27  which  represent  the  left,  central  and  right  parts 
of  the  model  respectively.  Flat  residual  moveouts  in  Figures  6.25,  6.27  and  the  upper 
part  of  Figure  6.26  indicate  the  correctness  of  the  estimated  velocity  in  these  areas.  The 
bottom  part  of  Figure  6.26  shows  incoherent  signals  that  affect  the  accuracy  of  velocity 
estimation  in  this  particular  area.  Notice  that,  at  the  same  common  image  gathers  (see 
Figure  6.28),  migrated  data  using  the  true  velocity  also  contain  obvious  incoherency 
in  residual  moveouts,  although  the  stacked  section  shows  a  good  image.  This  example 
demonstrates  that  for  an  extremely  complex  structure  it  is  very  difficult  to  identify  the 
correct  velocity  model  based  on  the  criterion  of  kinematic  coherence. 

6.4  Summary 

Imaging  complex  structures  (such  as  the  Marmousi  data)  requires  powerful  prestack 
depth  migration  algorithms  as  well  as  advanced  velocity  analysis  techniques.  The  per¬ 
turbation  method  in  this  thesis  provides  a  useful  tool  for  updating  a  velocity  model  by 
matching  a  criterion  based  on  prestack  information,  which  is  one  of  the  key  elements  in 
velocity  model  determination  (Versteeg,  1994).  In  order  to  obtain  a  more  satisfactory 
velocity  analysis  result  for  a  complicated  model,  as  Versteeg  concluded,  one  also  needs 
related  geologic  information  as  constraints. 
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Fig.  6.1.  The  true  velocity  model.  The  velocity  unit  is  m/s. 
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Fig.  6.2.  Synthetic  data:  (a)  with  offset  of  100  meters  and  (b)  with  offset  of  900 

meters. 
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PrestcLck  migration  with  the  initial  constant  velocity.  The  offset  is  100  meters 
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Fig.  6.4.  Two  common-image  gathers  from  the  migrated  data  with  the  initial  velocity. 
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Common  image  gathers 


Fig.  6.5.  Four  common-image  gathers  from  the  migrated  data  with  the  initial  velocity. 
The  numbers  1,  2,  3,  etc,  correspond  to  outputs  of  different  CIGs. 


Fig.  6.6.  Final  estimated  model  from  velocity  analysis. 
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Common  image  gathers 


Fig.  6.7.  Common  image  gathers  from  the  migrated  data  with  the  velocity  in  Figure 

5.6. 
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Fig.  6.8.  Migration  with  the  estimated  velocity  model  in  Figure  6.6. 
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Fig.  6.9.  Migration  with  true  velocity  model  in  Figure  6.1. 


Fig.  6.10.  Measured  velocity  model. 
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Fig.  6.11.  Physical-tank  data:  (a)  with  offset  of  268.2  meters  and  (b)  with  offset  of 

1243.6  meters. 
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cl  =  2438, 3962, 3309, 3570  x=1 280in  c2  s  3570, 51 82, 4779  x=1 280m 


Common  imajis  galhora  Common  image  gathers 


c3  =  4779, 7315, 6681  x=1280m  c4  s  6681 , 3962, 4869  x=3109m 


Fig.  6.12.  Velocity  analysis  on  the  physical-tank  data.  The  measured  velocities: 
Cj  =  3581,  C2  =  4801,  C3  =  6831,  C4  =  4801. 
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Fig.  6.13.  Estimated  velocity  model  from  velocity  analysis. 
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Fig.  6.14.  Five-olfset  stacked  migration  output  for  the  physical-tank  data,  AGC 
applied  to  the  stack.  The  input  velocity  model  is  shown  in  Figure  6.13. 
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Fig.  6.15.  Five-oflFset  stacked  inversion  output  for  the  physical-tank  data,  AGC  applied 
to  the  stack.  The  input  velocity  model  is  shown  in  Figure  6.1. 


0 


Midpoint  (km) 

Fig.  6.16.  Estimated  velocity  model  from  velocity  analysis  by  perturbation. 
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Fig.  6.17.  Five-offset  stacked  migration  output  for  the  physical-tank  data.  The  input 

velocity  model  is  in  Figure  6.16. 
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Fig.  6.18.  The  Marmousi  velocity  model.  The  darker  shading  denotes  higher  velocity 
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Fig.  6.19.  The  minimum-olFset  Marmousi  data.  The  offset  is  200  meters. 
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Fig.  6.21.  The  updated  velocity  model  using  the  paraxial  ray  tracer. 
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Fig.  6.22.  Comparison  of  the  true  velocity  model  to  the  estimated  one.  The  dark  curve 
denotes  the  true  velocity  and  the  gray  denotes  the  estimated  velocity.  The  top  figure  is 
at  location  x=8  km;  the  middle,  at  location  x=6  km;  the  bottom,  at  location  x=4  km. 
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Fig.  6.23.  19-ofFset  stacked  migration  output  for  the  Marmousi  data.  The  input 

velocity  is  one  in  Figure  6.21. 


Fig.  6.24.  19-ofFset  stacked  migration  output  for  the  Marmousi  data.  The  true  velocity 

is  used. 
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!5.  Ten  common  image  gathers.  19  offsets  in  each  CIG.  The  image  location 
ranges  from  4  km  to  4.25  km. 


Fig.  6.26.  Ten  common  image  gathers.  19  offsets  in  each  CIG.  The  image  location 

ranges  from  6  km  to  6.25  km. 
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17.  Ten  common  image  gathers.  19  offsets  in  each  CIG.  The  image  location 
ranges  from  8  km  to  8.25  km. 
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Fig.  6.28.  Ten  common  image  gathers  from  migrated  data  using  the  true  velocity.  The 
image  location  ranges  from  6  km  to  6.25  km. 
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Chapter  7 


CONCLUSIONS 


PrestcLck  depth  migration  provides  a  powerful  tool  for  knowing  if  a  migration  veloc¬ 
ity  is  acceptable,  because  of  its  high  sensitivity  to  a  velocity  model.  Moreover,  residual 
moveout  in  migrated  data  can  be  used  to  update  an  unacceptable  velocity.  In  this  thesis, 
I  studied  the  relationship  between  residual  moveout  and  migration  velocity  in  a  general 
context  and  derived  analytical  formulas  to  update  velocities.  For  handling  complex  me¬ 
dia,  my  formula  estimates  the  update  in  velocity  by  computing  a  derivative  function  of 
imaged  depths  with  respect  to  velocity  in  a  general  background  medium  context.  This 
formula  is  more  accurate  than  conventional  formulas  based  on  hyperbolic  residual  move¬ 
out  when  the  medium  has  strongly  lateral  velocity  variations.  In  this  sense,  I  made 
contributions  to  theory  and  applications  of  migration  velocity  analysis.  However,  some 
questions  remain  with  regard  to  efficient  application  of  migration  velocity  analysis. 

1.  Efficient  and  accurate  migration  algorithms 

In  order  to  handle  complex  structures,  repeated  prestack  migration  is  unavoidable 
in  velocity  updating.  Therefore,  the  efficiency  of  migration  velocity  analysis  largely  de¬ 
pends  on  which  prestack  migration  algorithm  is  used.  In  general,  integral-type  migration 
approaches,  such  as  Kirchhoff  or  Gaussian  beam,  are  preferable,  because  those  methods 
can  implement  either  common  shot  gathers  or  common  offset  gathers,  and  have  the  flex¬ 
ibility  to  image  the  targeted  structures  in  which  velocity  is  being  estimated.  Besides  the 
efficiency,  the  accuracy  should  be  considered  as  well.  For  imaging  complex  structures,  a 
migration  algorithm  should  be  designed  to  handle  turning  waves,  and  multivalued  trav- 
eltimes. 

In  the  Kirchhoff  integral,  traveltime  calculation  by  ray  tracing  or  flnite  differencing 
dominates  the  total  cost.  Finite  differencing  only  calculates  the  flrst  arrival  traveltime,  so 
this  approach  works  fast  but  fails  in  calculation  of  multivalued  traveltimes.  Ray  tracing 
by  a  shooting  method  may  handle  multivalued  traveltimes,  but  the  current  algorithm  is 
extremely  time  consuming  (Sun,  1992).  The  paraxial  ray  method  shows  advantages  in 
handling  multivalued  traveltimes  and  caustics.  This  method  uses  information  from  the 
standard  dynamic  ray-tracing  method  to  extrapolate  traveltimes  and  ray  amplitudes  at 
receivers  in  the  vicinity  of  a  central  ray.  In  general,  the  traveltime  calculation  by  the 
paraxial  ray  tracer  is  more  costly  than  by  flnite  differencing. 

The  cost  of  traveltime  calculation  in  Kirchhoff  migration  may  be  reduced  by  gen¬ 
erating  a  table  of  traveltimes,  if  the  computer  memory  and  disk  space  are  allowable.  In 
this  way,  the  traveltime  calculation  at  one  location  is  implemented  only  one  time,  so  the 
cost  of  traveltime  calculation  makes  up  relatively  small  portion  of  the  total  migration  cost. 
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2.  Interactive  computational  techniques 

In  practice,  a  velocity  estimation  procedure  must  involve  interactive  computational 
techniques,  which  is  especially  important  in  the  3-D  case.  During  velocity  analysis,  one 
needs  good  computer  graphics  tools  to  pick  imaged  depths  and  to  build  and  visualize 
the  velocity  model.  Because  of  the  limited  resolution  of  velocity  analysis,  mathematical 
estimation  of  velocity  models  may  be  unreasonable  and  non-unique.  By  means  of  com¬ 
puter  graphical  software,  one  can  interpret  and  modify  the  velocity  model  based  on  some 
geological  knowledge  and  information  from  other  sources. 
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Appendix  A 


DEPTH  VARIATION  TO  VELOCITY 


The  purpose  of  this  appendix  is  to  show  that  the  change  in  image  depth  z,  with 
respect  to  changes  in  velocity  is  positive  for  downward  propagating  waves  and  negative 
for  upward  propagating  turned  rays.  Suppose  that  migration  velocity  has  a  variation 
6c  =  6c(x).  For  a  fixed  horizontal  location  x,  the  imaged  depth  z  and  the  surface 
position  parameter  ^  are  functionals  of  migration  velocity.  One  takes  the  total  variation 
in  equation  (2.5)  with  respect  to  6c;  then 


‘dr^  dTr 


6^  +  [6cTg  +  6cTr 


^  dt 
6z  =  -H. 


(A.1) 


Here  6c  is  used  to  denote  the  partial  functional  variation  from  6c.  For  example,  6cTs 
means  the  variation  of  r,  related  to  changing  c  for  fixed  x  and  Using  equation  (2.6) 
to  elimate  dT^jd^  and  dr^jd^  yields 


dtr 
dz  dz 


6z  =  -[4r5(X5,  x)  +  6cTr{Xs,  «)]. 


The  derivatives  of  traveltime  with  respect  to  depth  can  be  represented  by 

9Ts{Xs,x)  _  C0S(6^)  dTr{Xr,x)  _  COs{6r) 

9z  c(x)  ’  dz  ~  c(x)  ’ 

where  6^  or  6r  are  the  angles  between  the  ray  direction  from  the  source  or  the  receiver, 
respectively,  and  the  downward  pointing  vertical  direction  at  x,  or,  equivalently,  the  angle 
of  incidence  from  the  source  or  receiver,  respectively,  measured  to  the  upward  vertical 
direction.  I  then  rewrite  6z  as 


6z  = 


c(x)[($cr3(x^,x)  +  6cTr{x^,x)\ 
COs{6i)  +  COs{9r) 


(A.2) 


An  increase  in  c—6c  positive— always  results  in  a  decrease  in  the  traveltimes— <5cr,+4T,.— 
negative,  so  that  the  numerator  on  the  right  is  negative.  More  precisely,  an  increase  in 
the  averaged  c  along  raypaths  does  the  same  work.  In  addition. 


COS(03)  +COs(0r)  =  2  cos 


78 


In  this  equation,  Or  —  Os  is  the  angle  between  the  two  raypaths,  so  \0r  —  ^a|/2  <  90°,  i.e.. 


cos 


'Or  -Os 


>  0. 


{0r+0s)l2  is  the  angle  between  the  normal  of  reflection  and  the  up- vertical.  For  downward 
propagating  waves,  {Or  +  Os) {2  <  90°,  so 


cos 


Or  +  Os 


>0; 


for  upward  propagating  waves  (turned  waves),  {Or  +  Os) /2  >  90°,  so 

^  Or  +  Os^ 


cos 


<  0. 


Thus,  when  taking  account  of  the  minus  sign  on  the  right  side  of  the  equation  (A. 2), 
I  conclude  that  the  variation  in  the  imaged  depth  with  respect  to  increment  in  the 
propagation  speed  is  positive  for  downward  propagating  waves  and  is  negative  for  upward 
propagating  (turned  waves).  This  is  true  for  any  velocity  distribution. 
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Appendix  B 


RMO  VELOCITY  REPRESENTATION  FOR  V(Z) 


In  this  appendix,  the  representation  of  RMO  velocity  is  derived  for  a  laterally  in¬ 
variant  velocity  medium  and  a  reflector  dip. 

From  Taylor’s  expansion, 


(l-p^v^(a))  ^^^  =  1  +  -v^(cr)p^  + 


so  that 


3 

/o .  2 

Using  equations  (B.l)  and  (B.2),  and  the  definition  in  (4.53),  I  obtain 


2  r*(nf -1-3/2  2,3.4  4.  2  .0/  4\ 

V —  =  — - -  -  o  ...  =v^  +  -(v;  -  vl)p^  4-  0{p^). 


r*(H- 3/2  n|p2) 

Thus,  the  derivation  of  equation  (4.55)  is  completed. 


(B.l) 


/  <i‘((!){l-p'‘v‘{a))  ^l'‘da=  j  v^da  =  T'(vl  +  ^vip'‘),  (B.2) 

^0  Jo  Z  Jo  I 

^  (l-pV(a))~^/2j^  =  r*(l-|-|t;|p2)  (g  3) 


(B.4) 
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Appendix  C 


LATERAL  VELOCITY  ANOMALY 


In  this  appendix,  the  representation  of  RMO  velocity  is  derived  for  horizontal  re¬ 
flectors  and  the  medium  with  a  lateral  velocity  anomaly. 

Suppose  that  reflectors  are  horizontal  and  the  true  slowness  w{x,z)  can  be  written 
by 

w{x,  z)  =  w{z)  (1  -I-  a{x,  z))  =  (1  +  a{x,  z)),  (C.l) 

V[Z) 

where  w{z)  is  a  reference  slowness  and  a{x,z)  is  a  small  perturbation.  The  two-way 
traveltime  t  in  the  medium  of  w  is  represented  by 

i{y,  h)  =  2  f  (C.2) 

where  9  is  the  angle  of  the  ray  path  from  the  vertical,  y  is  the  midpoint.  Because  the 
reflectors  are  horizontal  and  the  medium  velocity  is  a  depth-dependent  function  plus  a 
small  lateral  perturbation,  I  conclude  that  the  reflection  point  equals  approximately  the 
midpoint,  i.e,  x  =  y.  Under  the  assumption  of  a  small  perturbation,  the  ray  path  in  the 
medium  of  w{x,  z)  is  the  same  as  that  in  the  medium  of  w,  so  the  two-way  traveltime  in 
the  medium  of  w{x,  z)  is 


cos  9 


a{^{a),a)  da 
cos  9  v{a)' 


where  (^,  a)  is  a  point  on  the  ray  path  and 


(C.3) 


^(<t)  =  y  + 


tSiXi  9  ds. 


(C.4) 


In  particular,  ^(0)  =  h—y.  Taking  a  second-derivative  with  respect  to  h  in  equation 
gives 


d^t(h,  y) 

_  dH{h,y) 

CM 

+ 

o 

II 

da 

dh2 

1 

o 

II 

-c 

0  dh? 

cos  9 

h=o 

(C.3) 

(C.6) 


Furthermore, 


dh^  cos  9 


d‘^a  1  ^da  d  /  1  \  /  1  \ 

dK^  cos ^  j  ^  [^)  ' 


(C.6) 
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The  fact  that  0  =  0  at  /i  =  0  gives 


-f— )  = 

dh  \cos9j 


sin  0  _ 

cos^Odhj  ,  „ 


(C.7) 


COS0  |h=0=  1. 

Using  equations  (C.7)  and  (C.8),  I  simplify  equation  (C.6)  to 


dh^  cos  6 


d‘^a\  _  ,  ,5^/1 


h=o  /i=0 


Notice  that  ^  is  a  function  of  h,  so  it  follows  from  the  chain  rule  that 

d^a  _  d^a  /  V  9a 

df^  ""W  \^) 

From  (C.4),  ^  y  is  an  odd  function  of  h;  hence 

aft"  H.0  H.U  ■ 

This  result  and  equation  (C.IO)  give 

d^a  _  8“^ a  (9^^ 

9h‘^  ,=o~ 

Here  I  use  the  fact  that  ^  ^  at  h  =  0.  From  0  =  0  at  h  =  0,  and 

5^  /  1  \  /  sin 0  5^0 \  1  +  sin^ 0  ( d9\^ 


/  1  \  /  sm0  0^9 

Vcos0/  \cos^9dh‘^ 


cos^  0  \dh  )  ' 


I  obtain 


dh^  Vcos 


Vcos0/|^^jj  \dh  J 


In  order  to  do  further  calculations,  I  introduce  the  slope  parameter 


which  is  independent  of  a;  then 


(C.8) 


(C.9) 


(C.IO) 


(C.ll) 


(C.12) 


(C.I3) 


(C.14) 


(C.15) 


cosO  =  J1  —  (t;p)2. 


(C.16) 
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From  this  result  and  equation  (C.4), 

=  y  +  J  j  ds, 

V 1  “  {vpY 


and 


These  formulas  imply  that 


_  n  vp 
•/o  yjl  -  {vp)' 


■-ds. 


(C.17) 


(C.18) 


—  =  / 
dn  JtT 


dp  JtT  [1  — 


dp 


=  i 


[1  -  (up)2]3/2 

do  V 

dp  cos  0 ' 


ds, 


ds, 


Equations  (C.19)  and  (C.20)  yields 


'ae  fdhy^' 

dh 

h=0 

dp  \dp  J 

J  h=0 


yds  _  2  yds 
Jq  vds  ~  tov^{z)  ' 


(C.19) 

(C.20) 

(C.21) 

(C.22) 


where  is  the  asymptotic  stacking  velocity  of  v  and  to  =  ^(y,  0);  equations  (C.20) 
and  (C.21)  yield 


dO 

'do  fdhy^' 

dh 

/i=0 

dp  \dp  J 

v{a)  2v{a) 


J  h=0 


fovds  tov^iz)' 


(C.23) 


Substituting  equations  (C.12),  (C.14),  (C.22)  and  (C.23)  into  equation  (C.9),  I  obtain 


dh^  \  cos  0 


d^a  /2  vds\  .  4v^((7)  . 


Substituting  equation  (C.24)  into  equation  (C.7),  yields 


d^t 

dh? 


h=0 


dH 

dh? 


+ 


8  r  \d‘^a  {J^vds)^  av^{a) 


h=0 


r 

Jo  difi  V: 


+ 


^v‘l{z)  Jo  [dy"^  v‘l{z)  v‘^{z)  J  v{a) 


da 


(C.25) 


The  asymptotic  stacking  velocity  and  Vs  are  calculated  by 


w^{y,z)  = 


1 


^f(y,  2r) 


t{y,h) 


d'^t{y,  h) 
dh^ 


(C.26) 


h=0 
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Using  equations  (C.3),  (C.26),  (C.27),  and  notation  fo=^(2/,  0))  I  obtain 


(C.27) 


<(y,  0) 


dH{h,y)i 


-  t{y,  0) 


dH{h,  y) 


=  4^3  (y,  z)  - 

to 

.N  (^)  r  _  da 


=  4wf  (y,  z)  -  4w^iz)  -  ^  ^  a(y,  (C.28) 


that  is, 


>/-  \  -2/  \  r  t  \  «cr 

c{y,c)^  + 


da  t{y,  0)  r  dH 


{a)  4 


(C.29) 


When  a  is  small,  t  is  an  approximation  of  t.  By  ignoring  the  second-order  perturbation, 
I  obtain 


/,  0)  ’ 

4 


'  dH 

dH  ■ 

_  ^0 

'dH 

dH  ' 

».o 

» J 

~  4 

U'*"  fcO 

» J 

Using  equations  (C.25),  (C.29)  and  (C.30),  I  obtain 

,„2^,  _  2^  f\.dcr  ,  2  /•-  (/;  vds)^  ,  0:^2 (^r)]  rfa 


^i^^(y,2r)  -  w;(2)  = 


da  2  ’ 

Vo  v{a)  iov^  Jo 


ay2  t;2(2)  ^2(2)  J  v{a) 


(C.30) 


.  (C.31) 


Notice  that  y  =  x,  so  the  above  equation  becomes 

-,\  n  ..  r  o2  _  /  rz  —  I  \  2 


=  f  r  +a(.,.)(l  +  ^)  4^.,  (C.32) 

w^{z)  to  Jo  ydx^  \vs{z)  J  v^^{z)’^v{a)  ' 

which  is  the  result  I  stated  in  equation  (4.65). 
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Appendix  D 


EQUATION  SYSTEM  FOR  PERTURBATION 


In  this  appendix,  a  linear  equation  system  is  derived  for  parameter  estimation  in  the 
perturbation  method  of  Chapter  5. 

By  using  the  matrix  and  vector  notations  in  equation  (5.18),  I  obtain 


AiW  =  —  ■ 


"•  “i 


m  ^  r-r 

J  =  1  |=:1 


=  E 


=  E(aa,#), 


(D.l) 


zf’+Azf -iW+AzW  =  zf’-jW  +  Azf’-AiW 

=  #>-i5)+EAA,(s;f-i®).  (D.2) 


Define  a  quadratic  objective  function 


/(AA)  =  E  E  . 

fc=ii=i 

Thus,  finding  the  minimum  of  /(AA)  is  equivalent  to  setting  its  gradient  with  respect  to 
AA  equal  0,  i.e., 

5AA/  /-l,2,...,n.  (D.3) 


By  using  equation  (D.2), 


Km?  _  n  /  _ “I  * 

/(AA)  =  2:  zf  -  m + x;  AA,  -  #') 

k=lj=l  L  :=1  . 
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The  derivatives  computed  from  function  are 


3/(AA)  _  ^  ” 

SAA,  tih 


2 EE  [-f  - ^  +  E  AA,  {sir  -  #')]  (jJ?'  -  j<‘>) 

2EE(#’-5®)(fff-^ 

it=l  i=l  ^ 

+2  E  E  E  AA,  (slf  -  ^  (si!'  -  iPO  ■ 

Ar=:lj=:li=l  ^  ^  ^ 


=  2EE(#’-5®)K’-ft“’ 

it=l  i  =  l 


K  m  n 


(D.4) 


The  condition  (D.3)  implies  that 


K  m  n 


E  E  E  AA,  (s!‘>  -  ^  -  9|‘>)  =  -  E  E  (f  ’  -  (sff  -  #’) ; 

fc=lj=ii=l  /V  /  k=lj=l  ^  ^ 


E^'*' 

.ib=:l  « 


(D.5) 


Thus,  I  completes  the  verification  that  the  solution  of  equation  (5.21)  is  just  the 
minimum  of  the  left  side  in  (5.20). 
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Appendix  E 


DERIVATIVE  FUNCTIONS  IN  SIMPLE  CASES 


In  this  appendix,  the  representation  of  derivatives  of  imaged  depths  with  respect  to 
velocity  and  velocity  vertical  gradient  is  derived  for  a  horizontal  reflector  and  a  constant 
background  velocity. 

From  equation  (5.23), 


(E.l) 


where  L  is  the  raypath.  Here,  to  simplify  notation,  I  suppose  that  all  derivatives  with 
respect  to  A  are  evaluated  at  Ai  =  A2  =  0.  When  the  background  velocity  is  a  constant 
uo  as  in  equation  (5.24),  equation  (E.l)  becomes 


^  = _ L  r 

dX  cose  Jo  “  dX  ’ 


(E.2) 


where  9  is  the  angle  between  the  raypath  and  the  vertical.  For  the  velocity  function  in 
equation  (5.28),  if  .2:  >  d,  then 

dv{x,z')  _ 
dXi 


dv{x,  z') 


=  z  -  zq; 


if  z  <  d,  then 


dv{x,z')  dv{x,z') 


Therefore, 


or  ^ _ j_  r  ^  _ 

dAi  cos  9  Jd  Vq 

1  r  ^JZJldz'  = _ 

20s  9  Jd  Vq  cos  9 


1  z  —  d 
cos  9 


=  ~  ~  ■  (E.4) 

0X2  COS  9  Jd  Vq  cos  9  2Vq  '  ^ 

Because  the  reflector  is  horizontal,  the  ray  paths  from  the  source  and  the  receiver  are 
symmetric.  I  obtain 


cos  da  =  cos  dr  = 
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Substituting  the  above  formula  and  equations  (E.3)  and  (E.4)  into  equation  (5.13),  yields 

h?  +  z  —  d 


gi{x,h)  = 


and 


^  +  z‘^{z-  ZqY  -{d-  ZqY 

g2{x,h)  = - - - 


2vo 

Thus,  I  complete  the  verifications  of  equations  (5.29)  and  (5.30). 
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